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ABSTRACT OF DISSERTATION 


Heat Pipe Vapor Dynamics 
by 

Farrokh Issacci 

Doctor of Philosophy in Mechanical Engineering 
University of California, Los Angeles, 1990 
Professor Ivan Catton, Co-chair 
Professor Nasr M. Ghoniem, Co-chair 

The dynamic behavior of the vapor flow in heat pipes is 
investigated at startup and during operational transients. The 
vapor is modeled as two-dimensional, compressible viscous flow in 
an enclosure with inflow and outflow boundary conditions. 

For steady-state and operating transients, The SIMPLER 
method is used. In this method a control volume approach is 
employed on a staggered grid which makes the scheme very stable. 
It is shown that for relatively low input heat fluxes the 
compressibility of the vapor flow is low and the SIMPLER scheme is 
suitable for the study of transient vapor dynamics. 

When the input heat flux is high or the process under a 
startup operation starts at very low pressures and temperatures. 



the vapor is highly compressible and a shock wave is created in the 
evaporator. It is shown that for a wide range of input heat fluxes, 
the standard methods, including the SIMPLER scheme, are not 
suitable. A nonlinear filtering technique, along with the centered 
difference scheme, are then used for shock capturing as well as for 
the solution of the cell Reynolds-number problem. For high heat 
flux, the startup transient phase involves multiple shock reflections 
in the evaporator region. Each shock reflection causes a significant 
increase in the local pressure and a large pressure drop along the 
heat pipe. Furthermore, shock reflections cause flow reversal in the 
evaporation region and flow circulations in the adiabatic region. The 
maximum and maximum-averaged pressure drops in different 
sections of the heat pipe oscillate periodically with time because of 
multiple shock reflections. The pressure drop converges to a 
constant value at steady state. However, it is significantly higher 
than its steady-state value at the initiation of the startup transient. 

The time for the vapor core to reach steady-state condition 
depends on the input heat flux, the heat pipe geometry, the 
working fluid, and the condenser conditions. However, the vapor 
transient time, for an Na-filled heat pipe is on the order of seconds. 
Depending on the time constant for the overall system, the vapor 
transient time may be very short. Therefore, the vapor core may be 
assumed to be quasi-steady in the transient analysis of a heat pipe 
operation. 



Chapter I 


INTRODUCTION 


1.1 SCOPE 

A heat pipe is a device which transfers heat by evaporating 
and condensing a fluid. A wick structure is used to return the 
liquid to the evaporation zone by capillary effect. Figure 1.1 shows 
a typical sketch of a heat pipe which consists of the case, the wick 
structure and the vapor space. The input heat at one end of the 
pipe increases the temperature of the liquid in the wick structure. 
The liquid then evaporates, the pressure in the vapor core 
increases and the vapor flows to the other end. At this end the 
temperature is lower and the vapor condenses and releases heat. 
The liquid, then, flows from the condenser to the evaporator by 
capillary action through the wick structure. 

The heat pipe is considered as a heat transfer device with 
very high conductivity. It can operate in weightless environments 
and transfer a huge amount of heat under conditions of very low 
temperature differences. These advantages make the heat pipe 
unique in space applications, electronic cooling systems and fusion 
energy systems. Although the steady state operation of heat pipes 
is well established, the transient behavior is not entirely 
understood. 
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Fig. 1.1. Schematic of a heat pipe 

The reliability of a heat pipe for space applications depends 
on its ability to respond to a sudden change in the heat load, and 
startup and shut down of the thermal system. A comprehensive 
experimental and theoretical study is needed to understand and 
properly model the physical phenomena which occur in heat pipes 
during operational transients. 

1.2 HEAT PIPE APPLICATIONS 

Gaugler [1] was the first to develop the concept of a heat pipe 
in 1942 proposing to apply heat pipes for cooling the interior of an 
ice box. During the early sixties, intensive work in this field was 
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started by Grover and his coworkers [2] at the Los Alamos 
National Laboratory. Their objective was to apply heat pipes in 
spacecraft designs. The first quantitative analysis of heat pipes 
was performed by Cotter in 1965 [3]. Since then there have been 
many advances in the theory, design and practice of heat pipes. 

The fact that heat pipes are able to operate under very low 
gravity conditions and that they do not include any moving parts, 
has motivated the continued development of the heat pipe for 
space applications. In these applications, the temperature is very 
high ( 1 000-2000°C ) and liquid metals (such as sodium and 
lithium) are used as working fluids. 

A heat pipe was proposed by Grover [2] to supply heat to the 
emitters of thermionic electrical generators and to remove heat 

from the collectors of these devices. In space based nuclear power 
systems, heat pipes are considered as a means of transferring heat 
from the generating point to the radiators, providing a uniform 
temperature on the radiator surfaces [4]. 

Recently heat pipes have been proposed to be used in 

internal cooling of the leading edges of reentry vehicles and 

hypersonic aircrafts to reduce the peak temperature and alleviate 
the thermal gradients at the edge [5]. In this case, the working 
fluid (a liquid metal) which is initially in the solid state, is melted 
by the reentry heat load and the heat pipe operation is initiated. 

The use of high temperature heat pipes as a heat transfer device 



has also been considered in rocket engines to cool and make the 
nozzles isothermal [6]. 

In case of a loss of coolant accident in a fusion reactor, it is 
desirable to ensure that no first wall melting or evaporation 
occurs. Such melting or evaporation of constituent elements can 
lead to the dispersal of radioactive waste hazards. In contrast to 
fission systems, a fusion power core usually contains enough 
thermal capacity to distribute the after-heat efficiently. Since the 
major component of decay after-heat is within a thin zone close to 
the first wall structure, a design incorporating heat pipes would 
result in an effective and fast redistribution of the decay after- 
heat. A heat pipe with a proper working fluid may be chosen such 
that the mechanism of heat transport through the pipe will start 
above a threshold design temperature. This would make the heat 
pipe operative only when needed, that is in the case of a loss of 
coolant accident. 

1.3 THEORETICAL BACKGROUND 

The theory of heat pipe operation is briefly outlined in this 
section. Pressure drop and heat flux in a heat pipe are discussed 
for engineering calculations along with the limits on the heat pipe 
operation. The section concludes with a discussion of the transient 
operation of heat pipes. The detailed explanation of heat pipe 
theory may be found in Dunn and Reay [7]. 
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1.3.1 Pressure Drop 

In order for heat pipes to operate, the capillary driving force 
in the wick, Ap cap , must overcome the pressure drop in the liquid 
phase, Api , in the vapor phase, Ap v , and the gravitational head, if 
any, Ap g . That is 


APcap^APl+bPi+bP' 


d.l) 


vapor 


The pressure difference across the interface of liquid and 
is 


2(7 cos# 

Ap = 

r 


(1.2) 


where a is the surface tension, 6 is the contact angle, and r is the 
radius of curvature. 

Figure 1.2 illustrates the operational conditions of the vapor- 
liquid interface. The radius of curvature in the condenser section is 
very large, whereas evaporation causes the radius of curvature at 
the evaporator menisci to decrease. The difference between the 
radii of curvature in condenser and evaporator, provides the 
capillary force. Thus, the capillary driving force, using Equation 
1.2, is 


*Pcap = 2o 


r cos 0' cos 6 C ^ 


(1.3) 
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VAPOR 



Fig. 1.2. Operational level of interface between liquid and vapor 

where the subscripts e and c refer to evaporator and condenser, 
respectively, and r is the radius of the capillary pore. Usually r c is 
very large with respect to r e . This yields 

2<xcos0 

&Pcap = ; (1.4) 

The pressure drop in the liquid phase, Api, depends mainly 
on the structure of the wick and the properties of the liquid. For 
engineering calculations, the liquid phase is modeled as a liquid 
flow in a porous medium. Since the mass flow in the evaporator 
and condenser is not constant, an effective length, l e ff, is defined 
for the porous medium. The effective length, through which the 
mass flow is constant, is approximated as 



<* = /. + (/. H )/2 


(1.5) 
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where l a is the length of the adiabatic section which is between the 
evaporator and condenser sections. l e and l c are the lengths of the 
evaporator and condenser. Thus, Darcy's law is used to relate the 
pressure drop in the liquid phase to the wick and liquid properties. 


i.e.. 


Ap, 


PM. 


( 1 . 6 ) 


where p/ and p/ are viscosity and density of the liquid, m is the 
liquid flow rate, A w is the wick cross sectional area and K is the 
wick permeability given by 


K = 


2 er h 2 


where e is the wick porosity, r ^ is the hydraulic radius and ///?£/ is 
the friction factor-Reynolds number product for laminar duct flow. 
The values of K for different types of wick structure are tabulated 
in Chi [8]. 

The pressure drop in the vapor phase, Ap v , consists of the 
pressure drop in the adiabatic section, Ap va , and the pressure drop 
in the evaporator, Ap ve , and in the condenser, Ap vc . That is, 

Ap v = Ap„ + Ap„ + Ap K ( li7 ) 
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In order to calculate the pressure drop in the adiabatic section, the 
vapor flow in this region is assumed to be laminar Poiseuille flow 
in a circular tube [7]. Then the pressure drop Ap va , is given by 


Ap, 


w , 4 


( 1 . 8 ) 


where r v is the radius of the vapor core. Applying a one- 
dimensional analysis Cotter [3] calculated the pressure drop in the 
evaporator and condenser. The momentum equation was solved 
for the pressure gradient written as a perturbation expansion in 
powers of 1 /Re, where Re is the radial Reynolds number. The 
calculations resulted in 


Ap„ = 



(1.9a) 


Ap yc = 



(1.9b) 
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Upon substituting Equations (1.8) and (1.9) in Equation (1.7), we 
get 


Ap v = 


%Pjnl a 

W 4 


+ 0.074 


m 


PS* 


( 1 - 10 ) 
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However, Busse [9], by using a modified Poiseuille velocity profile 
in Navier-Stokes equations, has shown that Ap v can be expressed 
as 


Ap ¥ = 


w 4 


(l.H) 


where l e ff is defined in Equation (1.5). 

The hydrostatic head, Ap g , is written as 

Ap t = gp,Lsm<p (1.12) 


where <p is the inclination angle of the heat pipe and L is the total 
length of the heat pipe. 

Upon substituting Equations (1.6), (1.11) and (1.12) into 
Equation (1.1), the circulation rate, m is given by 


m = 


N? C ap-m L sin<j? 


Hi , 

PM. TTp/ 4 


(1.13) 


1.3.2 Heat Flux 

Neglecting the sensible heat which is small compared to the 
latent heat, the axial heat flux rate, <7, is given by 

i = (1.14) 
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where, m is the circulation rate of the working fluid and hf g is the 
latent heat of evaporation. 

From Equations (1.13) and (1.14) the heat transfer rate can 
be expressed as a function of the capillary driving force, the 
geometry and the properties of the working fluid; i.e., 


. _ (&P cap -gp,L sin (p)h fl 
q - r o' 


/ 


•ff 


PFK Wv 


* j 


(1.15) 


1 . 3.3 Limitations 

There are several limitations on the output heat transported 
through heat pipes, during the steady state or transient operations. 
These limitations are briefly described below [8]. 

1. Viscous limit: In long pipes and at low temperatures, the vapor 
pressure is low and the effect of viscous friction on the vapor flow 
may dominate inertial forces. In this situation the circulation of the 
working fluid is limited, which consequently, limits the heat 
transfer through the pipe. 

2. Sonic limit: At low vapor pressures, the velocity of the vapor at 
the exit of the evaporator may reach the speed of sound. Then the 
evaporator cannot respond to further decrease in the condenser 


pressure (i.e., the vapor flow is choked), which limits the vapor 
flow rate. 

3. Capillary limit: A capillary structure is able to provide 

circulation of a given fluid up to a certain limit. This limit depends 

on the permeability of the wick structure and the properties of the 
working fluid. 

4. Entrainment limit : The vapor flow exerts a shear force on the 

liquid in the wick which flows opposite the direction of the vapor 
flow. If the shear force exceeds the resistive surface tension of the 
liquid, the vapor flow entrains small liquid droplets (Kelvin- 

Helmholtz instabilities). The entrainment of liquid increases the 
fluid circulation but not the heat transfer through the pipe. If the 
capillary force cannot accommodate the increased flow, dryout of 
the wick in the evaporator may occur. 

5. Boiling limit: At high temperatures, nucleate boiling may take 

place which produces vapor bubbles in the liquid layer. The 

bubbles may block the wick pores and decrease the vapor flow. 

Furthermore, presence of bubbles decreases the conduction of heat 
through the liquid layer which limits the heat transfer from the 

heat pipe shell to the liquid which is by conduction only. 

The above limits are fully explained in several textbooks (for 
example, Chi [8]) which explain how they are to be determined. 

These limitations dictate the operational envelope for a given 
design of a heat pipe, as shown schematically in Fig. 1.3 [2]. The 


Axial Heat Flux 



Temperature 


Fig. 1.3. Heat pipe operating limits 

shape of the operational envelope changes for different working 
fluids and different wick materials. 



1.3.4 Transient Operation 

Transient operation of heat pipes may be considered in 
terms of three different regimes: startup, operating transients and 
shut down. The following is a brief description of these regimes. 

Startup 

Heat pipe startup is a transient process through which the 
heat pipe starts its operation from an static condition. The working 
fluid of a heat pipe may initially be very cold or frozen in the wick. 
Under this condition, the vapor pressure is very low and the vapor 
flow is free molecular. The heat transferred from the evaporator 
melts the frozen working fluid in the evaporator and evaporation 
takes place at the liquid-vapor interface. The vapor pressure is 
initially very low and the vapor flow is highly accelerated. Then 
the flow may be choked at the exit from the evaporation zone if its 
velocity reaches the sonic velocity. By addition of heat, the frozen 
working fluid is completely melted and continuum vapor flow 
prevails throughout the vapor core. 

Behavior of the vapor flow during the startup of a heat pipe 
from a frozen state may be considered in three distinct phases: (i) 
free molecular flow throughout the vapor core, (ii) establishment 
of continuum flow in the evaporation zone with a front, moving 
towards the condenser with possibility of choked flow at the 
evaporator exit and (Hi) continuum flow of vapor throughout the 


vapor core. 



Operating Transients 

Operating transients involve power changes during nominal 
operating conditions. There are several transient conditions of 
interest. A heat pipe may operate in a low power stand-by mode 
but be required to quickly respond to full power operation. A heat 
pipe may also be required to perform under load-following 
conditions, which depend on changing demands of the power 
conversion system. A nuclear reactor control system may also 
require changes in heat generation that would have to be handled 
by the heat pipe system. 

The physics of interest for these situations depend on the 
initial conditions, the rate of change and magnitude of input power 
and condenser conditions. A heat pipe operating in stand-by mode 
at low vapor pressure may encounter the sonic limit while 
powering up. Large temperature gradients and complex vapor 
dynamics would develop. After passing through the sonic-limited 
operation, or for a heat pipe initially operating at high vapor 
pressure, the capillary, entrainment and boiling limits are of 
concern. 

Shut Down 

Shut-down of a heat pipe refers to the transient process 
occurring when the power of the thermal system is shut-off. The 
input heat of the heat pipe is suddenly reduced to zero and the 
pipe cools down to a static condition. There are several problems 
encountered with this mode of transient operation. After the 
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power is shut off, the working fluid of the heat pipe cools down 
• and eventually solidifies. This phase change may affect the wick 

structure. The screen sheets of the wick are very delicate 
structures with fine meshes. Freezing of the working fluid in the 
I screen holes may damage the wick structure by tearing the screen. 

Also if for some reason, the working fluid freezes mainly inside the 
screen holes, it will cause a major reduction in the driving force of 
I the liquid flow. 

Another possible problem is the appearance of cracks or 
bubbles in the solidified mass of the working fluid. The presence of 
I cracks or bubbles can decrease the conduction of heat in the 

solidified mass during the startup, thereby the subsequent startup 
time increases. Cracks can also cause strain on the structure of the 
I heat pipe and the screen and fatigue may become a serious 

problem after many shut downs. A proper heat pipe design must 
take into account the problems involved in the shut down 
I transient mode, in order to have a reliable and long lasting heat 

pipe. 

I 

1.4 LITERATURE REVIEW 

The heat pipe studies were initiated by the experimental 

^ work of Grover et al. [2] and the quantitative analysis of Cotter [3]. 

Initial qualitative experiments were performed on heat pipes with 
water and sodium as working fluids. Cotter considered a right 
^ circular heat pipe with large length-to-diameter ratio. He 
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performed a simple one-dimensional analysis on the vapor flow to 
develop the general basic theory for quantitative engineering 
calculations of heat pipe operations. 

The vapor flow in heat pipes is a complicated problem due to 
high nonlinearities of the governing equations and the inflow and 
outflow boundaries in the evaporator and condenser. Different 
approaches have been used to simplify the problem. In most of the 
work done so far, the vapor flow is analyzed under steady state 
conditions, in one- [10-13] and two-dimensions [14-24]. 

Banskton and Smith [14,15] developed numerical solutions 
for incompressible, steady state, two-dimensional vapor flow in a 
heat pipe. Their results show flow reversal in the condenser for 
radial Reynolds numbers larger than two. Indication of reverse 
flow is important in the design of a heat pipe in calculating the 
shear forces acting on the wick structure. Flow reversal may also 
affect the entrainment limit. They also found that the pressure 
profiles across the heat pipe were two dimensional. Their study 
recommended the solution of the complete (compressible) two 
dimensional Navier-Stokes equations to provide an accurate 
prediction of the pressure drop along the pipe. 

McDonald et al. [16] analyzed a vapor flow in an enclosure in 
the presence of air as a noncondensible gas. The coupled governing 
equations for steady state two-dimensional flow were numerically 
solved to study the effects of air on the evaporation and 
condensation rates of the water heat pipe. They concluded that the 



two-dimensionality of the flow field had a very small effect on the 
concentration of the noncondensible gas and on the temperature 
distributions in the vapor core. 

However, Rohani and Tien [17,18] showed that in liquid 
metal gas-loaded heat pipes, a one-dimensional analysis is not 
adequate. This conclusion was explained by averaging the two- 
dimensional results over the cross section, and comparing these 
with those of a one-dimensional analysis [15]. They performed a 
numerical analysis of the steady state two-dimensional heat and 
mass transfer in the vapor gas region of a gas loaded heat pipe. 
Their results show that the vapor flow in the vicinity of the vapor- 
gas front is negligibly small and they concluded that heat and mass 
are transferred only by conduction and diffusion at the front. But, 
in an analysis done by Peterson and Tien [19], it was shown that 
the natural convection and radial diffusion have a significant effect 
on the noncondensible gas distribution in the vicinity of the vapor- 
gas front. 

Gas controlled heat pipes (GCHP) were recently investigated 
by Galaktionov and his colleagues from the Moscow Energy 
Institute [20,21]. In their analysis, they included convection and 
diffusion across the front under heat loads considerably smaller 
than the sonic limit of heat pipe operation. Galaktionov et al. [20] 
presented a mathematical model for the incompressible, steady 
state, two dimensional problem for the heat and mass transfer in 
the region of the vapor-gas front in a GCHP. In their analysis, using 
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boundary layer assumptions in the evaporator zone, they applied 
an integral method and approximated the axial velocity by a 
fourth-order polynomial. 

The same problem was analyzed by Galaktionov and 

Trukhanova [21], taking into account the compressibility of the 
vapor. By defining a stream function for the velocity field, the 
momentum equation was written in terms of the vorticity and 

stream function. The governing equations were solved numerically 
by the Gauss-Seidel method. 

They also carried out experimental studies of a planar gas 
controlled heat pipe. The results show that the variations of 
temperature in axial and radial directions are quite comparable, 

which clearly demonstrates the two dimensional character of the 
problem. The velocity profiles show the reverse flows in the region 
of the vapor-gas front. 

Reverse flow in the condenser section was also observed by 
Ooijen and Hoogindoorn [22] who performed experimental and 

numerical studies on heat pipe vapor flow. Their analysis was 
confined to a steady-state, incompressible, two-dimensional flow. 
They also found that the total pressure drop over the heat pipe, for 
high radial Reynolds numbers, was higher than that approximated 
by a Poiseuille flow model. 

Narayana [23] performed a numerical study of 

incompressible, steady state, two-dimensional vapor flow. He 
showed that the pressure terms in the Navier-Stokes equations 
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could be eliminated by integrating the momentum equation in the 
axial direction, on the pipe cross section. In this analysis, using 
boundary layer assumptions, the radial pressure variation was 
neglected. 

A double-walled concentric heat pipe was proposed by 
Faghri [24]. Two pipes with different diameters create an annular 
space for the vapor flow. The wick structures are attached to the 
inner surface of the outer pipe and the outer surface of the inner 
pipe. In this kind of heat pipe the area of heat transfer into and 
out of the pipe is increased and its efficiency of performance is 
expected to be higher. Faghri has carried out a numerical study of 
the vapor flow in such a double-walled heat pipe. He used an 
implicit marching finite difference method to solve the 
incompressible, steady state, two dimensional flow problem. His 
results show that for low Reynolds numbers, viscous effects 
dominate. For radial Reynolds numbers greater than one, pressure 
decreases in the vapor flow direction in the evaporator section, 
whereas it increases in the condenser section due to partial 
dynamic pressure recovery from the decelerating flow. For high 
radial Reynolds numbers, flow reversal was shown to occur in the 
condenser section. 

In studies of the dynamic behavior of heat pipes, three 
different areas of concern are startup, shut-down and operational 
transients. Among these the startup mode of transient is the most 
difficult one. As mentioned above, the startup of a heat pipe may 


be encountered with melting of the working fluid at solid state, 
free molecular flow of vapor and dry out and rewetting of the 
wick. 

Extensive studies of startup and shut down operations of 
heat pipes [25*29] have been carried out at the Los Alamos 
National Laboratory. The objectives were the applications of high 
temperature, liquid metal heat pipes to space power system heat 

rejection. Experiments were performed on a long heat pipe with a 
molybdenum container. Using liquid sodium and potassium as 
working fluids, design characteristics of high temperature heat 
pipes have been evaluated. For the temperature range of 700 to 
1600°AT, the useable length of 30m with length-to-diameter ratio 
of 300 was found for potassium heat pipes, whereas for sodium 
heat pipes useful lengths of more than 40m with length-to- 

diameter ratio of 800 were predicted. 

In the numerical studies by Castello et al. [29], heat pipe 

performance was modeled to include startup from the frozen state 
of the working fluid, with provision for free molecule flow of the 
vapor. Vapor core dynamics were simulated by a transient one 
dimensional compressible flow with friction and mass flow from 

the wick. The governing equations were solved using the KACHINA 
method [30], a numerical code is still under development. In the 
present status it requires 30 minutes of the Cray computational 
time per hour of actual time. 


Bystrov and Goncharov [31] performed experimental and 
theoretical studies on a heat pipe with and without a foreign gas. 
The steady state operation and startup from a frozen state of 
liquid were analyzed. By averaging over the axial cross section, 
they solved unsteady, one dimensional governing equations for the 
vapor and liquid phases. Their numerical scheme could be run only 
for very small time steps (10‘ 5 to 10‘ 7 s) because of numerical 
stability problems and as a result needed a large amount of 
computational time. In order to increase the time steps, they used 
a lumped parameter approach for heat balance and avoided 
calculations of detailed temperature distributions. 

Figure 1.4 shows the qualitative features of the temperature 
fields for startup of a heat pipe with a noncondensible gas, 
predicted by Bystrov and Goncharov [31]. Sharp discontinuities in 
the temperature distribution are the result of shock waves at the 
deceleration of vapor from supersonic flow. The heat pipe being 
initially at a uniform temperature (curve 1), is heated up 
uniformly and melting takes place (curve 2). The temperature 
increases and an intensive evaporation begins. As the vapor flows 
into the condenser the velocity increases up to the local sonic 
speed (curves 3 and 4). The volume of the condenser occupied by 
vapor increases and the noncondensible gas front is compressed 
further into the condenser (curves 5, 6 and 7). Supersonic flow 
occurs over the length of the active condenser. Curves eight 
through ten show the shock reflection back to the evaporator. 




Fig. 1.4. Frontal startup temperature history of a gas-loaded 

heat pipe [31] 

Further heating raises the vapor pressure and eliminates the 
shock wave (curve 11). Finally, the heat pipe reaches near 
isothermal conditions at its full operating temperature (curves 12 
and 13). 

Experimental studies of the heat pipe startup were 
performed by Ivanovskii et al. [32]. They measured the 
temperature distribution along a sodium heat pipe initially at 
frozen state at ambient temperature. Temperature measurements 
were made in the vapor by a moveable microthermocouple placed 
directly in the vapor core. Free molecule flow and transition to 
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continuum flow were observed at the early stages of the 
experiment. Then the vapor flow reached the sonic limit and a 
shock wave was created which caused an abrupt temperature 
drop. Non-isothermal conditions during sonic limit operation were 
seen to exceed 300 °C . Continuum vapor flow was found to occupy 
the whole vapor core and vapor pressure increased along the pipe. 
The process then proceeded to a steady state condition with an 
isothermal temperature distribution along the pipe. The transient 
time from a frozen state to a steady state condition was reported 
to be about ninety minutes. 

Chang and Colwel [33-35] studied the transient operation of 
low temperature heat pipes experimentally and numerically. 
Experiments were performed with Freon 11 as a working fluid to 
study the dryout and rewetting of the liquid phase. The 
temperature range in experiments was between the room 
temperature and well above the critical temperature. At low heat 
fluxes, the liquid temperature was nearly uniform in the slab wick 
and only small gradient existed in the vapor. The measured 
transient time to reach steady state was about ninety minutes. At 
high heat fluxes, partial and complete drying of the wick were 
observed. In their analysis the vapor core temperature was 
assumed to vary with time only, allowing the solution of a two 
dimensional conduction problem for the liquid and wick structure. 

Ambrose [36] studied the pulsed startup of heat pipes. A 
pulsed startup refers to applying a step increase in input heat flux 
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to a heat pipe under steady state conditions. They performed 
experiments on a simple-screen-wicked copper-water heat pipe to 
investigate the transient time and the dryout and rewetting 
behavior. The transient time was found to be dependent on the 
capillary force, the cooling mechanism and the heat capacity of the 
pipe. Dryout of the wick at the evaporator was found to occur only 
for heat fluxes greater than the capillary limit. Dryout and 
rewetting were modeled using a lumped parameter method. Good 
agreement between the experimental data and the theoretical 
results was reported. 

Bowman [37] also studied the dynamics of vapor flow in heat 
pipes analytically and experimentally. In the experiments an air 
flow was pumped through a porous medium into a pipe at one end 
and sucked from the other end. The pressure and velocity field in 
the vapor core were measured. He also solved the steady state 
compressible, axisymmetric two-dimensional form of the Navier- 
Stokes equations numerically. Using the experimental and 
analytical results, he derived a friction coefficient expression to be 
used in a one dimensional transient model of heat pipe vapor flow. 
The friction coefficient is expressed as a function of local axial 
Reynolds number, Mach number, pipe aspect ratio and radial 
Reynolds number. 

Recently Jang [38] studied an overall analysis of heat pipe 
startup from the frozen state, numerically. In this analysis, the 
conjugate problem of conduction in the heat pipe case and melting 
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and convection in the liquid phase is solved during the transient 
mode. In the vapor core, the free molecule flow, choked and 
continuum flow of vapor are taken into account, separately. 
However, the vapor core is assumed to be at steady state at any 
flow regime. This is rationalized by the fact that the transient time 
for vapor flow to reach a steady state is much shorter than that for 
the liquid flow. 

1.5 OBJECTIVES 

The research performed so far on the transient behavior of 
heat pipes mainly concerns the startup and shut down operations. 
However, these modes of transients are still not well established. 
Furthermore, no major work has been done on the operational 
transient of heat pipes. For a comprehensive study of the dynamics 
of heat pipes, the vapor flow should be thoroughly investigated 
theoretically and experimentally. 

The previous theoretical studies of dynamics of vapor flow 
have considered a transient one-dimensional model. In these 
studies the friction coefficients of the vapor flow boundaries have 
been approximated using a steady state two-dimensional model. In 
addition, in steady state two-dimensional analyses it has been 
shown that the one-dimensional model of vapor flow is not able to 
accurately predict the axial heat and mass transfer and the axial 
pressure drops. Furthermore, two-dimensional steady-state 
studies indicate that flow reversal takes place in the condensation 



section. It is therefore important to establish the conditions for this 
mode during transient operations. 

The study of dynamic behavior of heat pipes at UCLA started 
by analyzing the liquid and vapor phases separately [39,40]. The 
objective was to develop very detailed computer codes for each 
phase, and then to combine them in a comprehensive code. Roche 
[39] developed a numerical code for the liquid phase during the 
startup transient mode by assuming the vapor phase at a constant 
temperature and pressure. He used a kinetic theory approach to 
model the phase change at the liquid-vapor interface. The model 
results in an expression which consists of a large coefficient 
multiplied by a very small temperature difference and requires 
very small time steps for stability. The large coefficient also 
transforms temperature differences that are otherwise beyond 
machine accuracy into computationally significant terms. He 
recommended investigation of an alternative model for the phase 
change process. 

The objectives of this work are to investigate thoroughly the 
dynamic behavior of the vapor flow in heat pipes at startup and 
during operating transients. The vapor flow in a heat pipe is 
mathematically modeled and the equations governing the 
transient behavior of the vapor core are solved numerically. The 
modeled vapor flow is a transient, two-dimensional, compressible, 
viscous flow in a closed chamber. The boundary conditions include 
inflow and outflow of heat and mass which simulate the 



boundaries of the vapor core in a heat pipe. The different solution 
methods are discussed and compared in Chapter Two. The 
numerical results of these methods are presented in Chapters 
Three through Five. 

The need for and difficulties of experiments on the vapor 
core are then discussed in Chapter Six. In this Chapter the proper 
optical methods for temperature and velocity measurements are 
described along with the advantages and shortcomings of these 
methods. This Chapter also includes a rather simple but 
fundamental experiment to study the dryout and rewet of the 
liquid phase. Conclusions and recommendations of this work are 
given in Chapter 7. 



Chapter II 

THEORETICAL MODEL OF THE VAPOR CORE 


2.1 INTRODUCTION 

Flow patterns found in a heat pipe are shown schematically 
in Fig. 2.1. The working fluid evaporates in the evaporation zone. 
The local pressure increases and the vapor flows towards the 
condensation zone, at which it condenses back to the liquid phase. 
The liquid flows from the condenser to the evaporator by capillary 
force through the wick structure. The main concern of this analysis 
is the vapor core response to changes in the evaporation and 
condensation rates due to a sudden increase, or decrease, in the 
input heat flux, or the condenser temperature. 

In this chapter the dynamic behavior of the vapor flow is 
analyzed using a transient two dimensional compressible flow 
model. The proper equations which govern the model problem will 
be written along with initial and boundary conditions which 

properly simulate heat pipe transient operations. Different 
numerical schemes and their advantages and disadvantages will 
then be discussed for different operational regimes. The 



I 


I 


I 


I 


I 


I 




v> 


Vapor Core 


Liquid Flow 





17 / 


Porous Wick 





A 


! t 

Input Heat Flux 
h — Evaporator — 


■Pipe Wall 


■Adiabatic Section 


Output Heat Flux 
- — Condenser — 1 


Fig. 2.1. Flow patterns in a heat pipe 

I 

application and numerical results of these schemes will be 
presented in Chapters Three through Five. 

I 

2.2 GOVERNING EQUATIONS 

Vapor flow in the heat pipe vapor core, shown in Fig. 2.1, is 
I modeled as a channel flow shown in Fig. 2.2. The bottom boundary 

of the channel is a thin porous medium which includes the liquid. 
The input heat flux to the evaporator and the temperature of the 


I 







• Fig. 2.2. Vapor core model of a heat pipe 

outer surface of the condenser are specified. The planar side walls 
I are assumed adiabatic. 

The equations governing the vapor flow are time dependent, 
viscous, compressible momentum, continuity and energy 
1 equations. An equation of state is used to relate pressure, density 

and temperature within the vapor core. These equations in 
Cartesian coordinates ( x’,y ') are 
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(2.1) 


( 2 . 2 ) 


( 2 . 3 ) 


( 2 . 4 ) 


( 2 . 5 ) 


In these equations, p' is density, u' and v' are velocities in x' and 

y' directions, respectively, T is temperature and p' is pressure. 
The properties p\ c' and k’ are dynamic viscosity, specific heat 

and thermal conductivity of vapor, respectively. 

In the momentum equations, 2.2 and 2.3, the last term in 
brackets accounts for compressibility of the fluid. Viscous 
dissipation in flows at high velocity or with high velocity gradients 
is not negligible. In heat pipe transient operations, where the 
vapor velocity may reach the speed of sound (sonic limit), with a 
large gradient at the initial stage of the process, the viscous 
dissipation is significant. This is taken care of by the last term in 
brackets in the energy equation 2.4. In the right hand side of the 
same equation, the terms involving pressure p', are presented 
owing to the compressibility of the flow. 

2.2.1 Normalization of the Equations 

The governing equations are normalized by using 
characteristic parameters of the problem, for parametric studies 
presented later. The characteristic length is the vapor core 
thickness, b. The characteristic velocity is related to the inflow 
mass flux of vapor in the evaporator, m"; that is 


where q " is the maximum input heat flux and hj g is the latent 

heat of vaporization. Then the characteristic velocity is 



(2.7) 


In this equation p' is density of the vapor at a reference condition 
discussed below. Having the characteristic length and velocity, 
time is normalized by t c , defined as 


t = A = £> 

tc u. 


m. 


(2.8) 


The normalized parameters are therefore defined as 


t x,y) = (x',y')/b , 
p = p / {p' r U c \ 


o u,v)=(u'y)/u e , 

T = RT'/U C 2 


t = t'U. / b. 


(2.9) 


and 

<p = (p'l (p' r , (2.10) 

where R is the gas constant. Here <p' is p\ p\ k', or c' and the 
subscript r stands for the property value at a reference condition 
according to the initial conditions of the process. 

The governing equations, 2.1 through 2.5, in nondimensional 
form, using normalized parameters, 2.9 and 2.10, are written as 
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P = f(pJ) 


(2.15) 


The nondimensional parameters are Reynolds number, Fe, Prandtl 
number, Pr, and specific heat ratio, y, defined as 





(2.16) 
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2.3 INITIAL AND BOUNDARY CONDITIONS 
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The working liquid is assumed initially at a very low 
temperature (close to the freezing point). Further there is no input 
heat and the stagnant vapor is in thermodynamic equilibrium with 
the liquid. The boundaries of the vapor core are shown in Fig. 2.2. 
A no-slip condition for the velocity and an adiabatic condition for 
the temperature are assumed on the side walls: 

@ x = 0 and x = L , 

u= 0 , v = 0 , — = 0. (2.17) 

dx 

On the center line the symmetry condition implies, 

@ y = b , 

X = 0 ,v =0 , f-o. (2.18) 

dy dy 

The boundary conditions at the liquid-vapor interface are 
the challenging ones. The liquid flow is assumed to be in a porous 
medium of thickness S which is much smaller than the vapor core 
thickness b. The axial velocity is assumed zero on this boundary. 
That is, 

@ y = 0 , u = 0 . (2.19) 
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To assign boundary conditions for the temperature and vertical 
velocity, the liquid-vapor interface is divided into three regions. 
In the evaporation zone the input flux, q", is a given parameter 
and the input flow is approximated by, 

@ y = 0 and 0 < x < L e , 

P v = ' h "~-r-^ ’ T =Tsat(P )■ (2.20) 

V T) 


In the above equation, hf g is the heat of vaporization and 
conduction in the liquid layer is neglected. The temperature is 
assumed to be the saturation temperature of the liquid 
corresponding to the interface pressure. In the adiabatic zone the 
boundary conditions are 
@ y = 0 and L e < x <{L e + L c ) , 

v = 0 , |^ = 0. (2.21) 

dy 

In the condensation zone the temperature of the outer 
surface, Tc is given. By equating the heat of evaporation to the 
heat conduction in the liquid layer, the outflow from this region is 
approximated by, 

@ y = 0 and ( L e + L c ) <x <L , 


pv = m" 



T-T 
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, T = T S at(P ) , 


( 2 . 22 ) 


where k e jj- is the effective conductivity of the liquid layer and the 

temperature is assumed to be the saturation temperature 
corresponding to the pressure at the interface. This neglects 
transient conduction in the fluid layer which may be an important 
contribution to the process. 

2.4 SOLUTION METHODS 

The five governing equations, (2.11) to (2.15) are solved 
numerically for the five variables, p, pu, pv, T and p. The boundary 
conditions associated with these equations are Eqs. (2.17) to (2.22). 
The boundary conditions on the vapor-liquid interface depend on 
the interface pressure, which is not known a priori. However, in 
the numerical methods described below, the pressure does not 
need to be specified on the boundaries, and is not calculated there. 
After completion of calculations, an interpolation method is used to 
determine the pressure on the boundaries. 

For low compressible flow, Issacci et al. [41] showed that the 
staggered grid made the numerical scheme very stable. In heat 
pipe transient operation, there are regions of high as well as low 
compressibility. Therefore, the governing equations are finite 
differenced on a staggered grid, shown in Fig. 2.3. The variables on 
this grid are not evaluated at the same point. By choosing control 
volumes which have the main grid points at their center, the 
velocity components are evaluated at the points that lie on the 



I Fig. 2.3 Staggered grid 

faces of the control volumes. The locations for u and v are shown 
I by short arrows in the figure. The other dependent variables, 

namely, p, T and p, are calculated for the cell centers shown by 
small circles. The cell centers do not lie on the boundaries. 
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Several different methods are available to solve the 
governing equations. Among these, the SIMPLER scheme, 
introduced by Patankar [42], is the simplest one. In this method a 
control volume approach is employed on a staggered grid which 
makes the scheme very stable. The main advantageous point is 
that the pressure field is guessed and the momentum equations 
are solved for the velocity field. Then the continuity equation is 
used to correct the pressure for the next iteration. Therefore for 
the problems where pressure is not specified on the boundaries, 
this scheme is very useful. When the continuity equation is used 
for pressure correction, the density is assumed to be independent 
of pressure. This limits the scheme to be limited to incompressible 
flows and low compressible flows. In Chapter Three the algorithm 
of this method is discussed in detail. Application of the SIMPLER 
scheme in the operating transient mode of the heat pipe vapor 
dynamics, is also presented in Chapter Three. It will be shown that 
for relatively low input heat fluxes the compressibility of the 
vapor flow is low and the SIMPLER scheme results in some 
interesting calculations. 

When the input heat flux is high or the process under a 
startup operation starts at a very low pressure and temperature, 
the vapor is highly compressible and a shock is created in the 
evaporator. Therefore, in the general analysis of heat pipes with a 
wide range of input heat flux, the standard methods, including the 
SIMPLER scheme, are not suitable. A nonlinear filtering technique, 
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along with the centered difference scheme, was used by Issacci et 
I al. [43] for shock capturing as well as the cell-Reynolds-number 

problem. The basic principles of this filter and its application in a 
one-dimensional model of heat pipe vapor dynamics will be the 
I subject of Chapter Four. It will be shown that the nonlinear 

filtering technique is quite simple to implement, is significantly 
more efficient than the standard methods, and is also found 
I suitable for calculation of shock reflections. Then the full 

governing equations of the two-dimensional model will be solved 
using the nonlinear filtering technique. The results of these 
I calculations are presented in Chapter Five for low- and high-input 

heat fluxes. 
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CHAPTER III 

OPERATING TRANSIENT MODE OF HEAT PIPES 
3.1 INTRODUCTION 

In thermal systems in which a heat pipe is used as a device 
to conduct a large amount of heat in a short time, the question 
"How much, how fast?" is the main issue in designing such systems. 
The research performed so far on heat pipe transient behavior, is 
mainly concerned with the startup and shut down operations and 
the analyses are carried out mostly on a one-dimensional model. 

In this chapter, the transient behavior of vapor flow is 
analyzed numerically using the two-dimensional model, described 
in Chapter Two. The main concern is the vapor core response to 
changes in the evaporation and condensation rates due to a sudden 
increase, or decrease, in the input heat flux, or the condensation 
temperature. 

In operating transient operation, when a step change in the 
input heat flux is imposed on a steady-state operation at moderate 
pressure and temperature, the vapor compressibility is low. Under 
these conditions the SIMPLER method is quite efficient for 
analyzing vapor flow. An algorithm for this method and its 
application in heat pipe transients are described below. The results 
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of calculations for a water heat pipe following a sudden change in 
input heat flux are then presented. 

3.2 THE SIMPLER METHOD 

The governing equations, (2.11) to (2.15), derived in Chapter 
Two, are solved for p, pu , pv, T and p. The boundary conditions are 
given by Eqs. (2.17) to (2.22). The equations are finite differenced 
on a staggered grid shown in Fig. 2.3. The method used here is 
based on the SIMPLER method described by Patankar [42]. The 
time-dependent terms are discretized using the backward Euler 
method while the power law scheme was used for differencing the 
convective terms. 

The calculational algorithm is as follows: The system 
pressure is assumed to be composed of two parts 

p{x,y,t) = p(t) + p(x,y,t) (3.i) 

where pit) is the space averaged system pressure whereas p(x,y,t) 
includes the space variation of the pressure. Assuming a pressure 
distribution, the momentum equations (Eqs. (2.12) and (2.13)) are 
solved for the velocity components u and v at a time step. Then 
using the continuity equation and density p from the last iteration 
on temperature, the pressure is corrected to be used again in the 
momentum equations. Upon accomplishing these iterations, the 
velocity field and p(x,y,t) will be found. 
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Using the velocity field, the energy equation (Eq. (2.14)) is 
• then solved for temperature. In order to calculate the density from 

the equation of state, and to evaluate the boundary conditions, we 
need the absolute value of the pressure. Global mass balance is 
I used to calculate the absolute pressure. The total mass, for an ideal 

gas, m , is related to the pressure by 

I m = jpdV = \-£-dV (3.2) 

V V 

where V is the volume. The increase in the total mass in one time 

step is 

I 

m(t + At) - m(t) = (m„ - rh oia )At (3.3) 


I 


I 


I 


I 


where m u and m oul are the input and output mass fluxes, 
respectively, and are calculated from the boundary conditions. 
Equations (3.2) and (3.3) are combined to yield 



f 


p\t+ts.t 



RT\ t ) 


dV 


(3.4) 


Thus Eq. (3.4) is used to evaluate the absolute pressure of the 
system. For a real gas, using an equation of state, the above 
procedure is followed to calculate the total pressure. 
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3.3 RESULTS 

The computational procedure described above was 
implemented in a two dimensional time-dependent compressible 
flow with the same type of heat pipe boundary conditions 
mentioned above. The computations were performed for several 
evaporator heat fluxes for water as a working fluid. The pipe 
dimensions used in the calculations were b = 2.5 cm and L = 3b 
and the computational grid spacings were Ax =0.1 and Ay = 0.05. 
The temperature and pressure variations, and consequently, the 
density variation, were found to be quite small inside the pipe. 
This is expected of systems in which condensation takes place 
without a non-condensible gas. 

The steady state results are shown in Figs. 3.1, 3.2 and 3.3. 
The steady state was assumed to be reached when the relative 
change of all variables in one time step was less than 0.5%. Figures 
3.1 show the axial velocity profiles at different axial cross sections 
along the pipe, for low-input heat flux, Fig. 3.1(a), and for high- 
input heat flux. Fig. 3.1(b). The flow stream lines for these two 
cases are shown schematically in Fig. 3.2. 

When the input heat flux is relatively low (/?e = 1.0), the 
evaporation rate is low and the vapor flows smoothly towards the 
condensation region (Figs. 3.1(a) and 3.2(a)). However, with high 
input heat flux (fle=1000) the evaporation rate is high, the vapor is 
ejected with a high momentum from the evaporator and under 
steady state conditions, condenses at a high rate at the condenser. 



Fig. 3.1 Axial velocity 
(b) high-i 



ofile. (a) Low-input heat flux (/?e=l), 
it heat flux (/?e=1000) 
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Qout 


(b) 


Fig. 3.2 Vapor flow pattern, (a) Low-input heat flux, 
(b) high-input heat flux 
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The high momentum flow at the two ends of the adiabatic 
section results in circulation flow in this region. That is, on the 
adiabatic surface there is a reverse flow. In addition, in the 
condensation region the high momentum vapor flow impacts the 
rigid end wall, returns and then condenses in the liquid region. 
This return flow is shown in Figs. 3.1(b) and 3.2(b). The vertical 
velocity profiles for different vertical cross sections are shown in 
Fig 3.3. The flow is clearly two dimensional. 



Fig. 3.3 Vertical velocity profiles 


It is worth noting that Galaktionov and Trukhanova [21] 

observed reverse flow in the gas-controlled heat pipe. To obtain 

the high heat flux solution, we used an incremental procedure 
where the magnitude of the input heat flux was increased and 
previous solutions were used as initial guesses for higher heat 
fluxes. It has been shown by Serrin [44] that the Navier-Stokes 
equations have a unique solution if they are considered as an 

initial value problem. Therefore, reverse flow is both unique and a 
valid solution for a high input heat flux. 

Transient behavior of vapor flow was studied by using a 
steady-state solution as the initial conditions for different input 
parameters. As an example. Fig. 3.4 shows the axial velocity 
profiles at the middle of the pipe, at different time steps. Curve 0- 
on this Figure is the steady state solution for Re = 100 where 

Reynolds number is defined as 


Re 


™,n b 


(3.5) 


and m u is the input mass flux. Curves 1 to 4 show this velocity 
profile at different time steps after a sudden increase in input heat 
flux for Re = 200. At the outset, the increase in input heat flux 
overcomes the reverse flow, see curve 1. As the velocity profile 
develops with time, the steady state profile, curve 4, shows a 
higher reverse flow at the liquid-vapor interface. 
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Fig. 3.4. Axial velocity profile at different time steps. 0- initial 
profile (t=0), 1- t=10* 4 , 2- t=2.5x10‘ 3 , 3- t=10‘ 2 , 

4- steady state (t=2.5x10‘ 2 ) where r=ta/b 2 . 

3.4 CONCLUSIONS 

The dynamic behavior of vapor flow in heat pipes was 
studied in this chapter for transient operations. Although the 
solution is presented for Cartesian coordinates, the computational 
code can be easily modified for circular heat pipes. Interesting 
phenomena have been observed in the adiabatic and condensation 
regions, in these sections high input heat flux causes flow reversal. 
This can result in negative shear forces on the wick structure. 
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Further, in transient operations the shear force at one location may 

• change sign. Detailed information on the magnitude and direction 

of the shear force is needed for design purposes, and this needs a 
proper two-dimensional analysis. 

• 
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Chapter IV 

NONLINEAR FILTERING TECHNIQUE 


4.1 INTRODUCTION 

Shock capturing and treatment of cell-Reynolds-number 
problems are among the major challenges in numerical 
approximations of nonlinear equations associated with fluid flow 
and heat and mass transfer. The earlier schemes proposed for 
shock computation, based on centered differencing, typically use 
artificial dissipation to control oscillations in the shock 
neighborhood. The only complete remedy is the use of monotone 
schemes, such as that by Godunov [45], which are only first order. 
The recent methods for shock capturing, using field-by-field 
decomposition and flux limiting (e.g., Harten [46] or Chakravarthy 
and Osher [47]), are able to capture a sharp discontinuity without 
oscillations. However, these methods are quite complicated and not 
efficient numerically. Filters have also been used in shock 
capturing [48], however a major problem is low accuracy and 
smearing of the shock front. Recently, Engquist, Lotstedt and 
Sjo green [49] proposed a nonlinear filter to accurately capture a 
shock without oscillations. The algorithms for this filtering 



technique are quite simple to implement, and they are significantly 
more efficient than earlier methods with comparable ability. 

The cell-Reynolds-number (or cell-Peclet-number) problem 
refers to a class of fluid dynamic problems with large-grid 
Reynolds number. Discretization of advective terms in these cases 
plays a significant role. The centered-difference scheme (CDS) is 
ill-behaved for grid Reynolds numbers larger than two because the 
coefficients of the difference equation are not necessarily positive, 
consequently, diagonal dominance of the corresponding system 
matrix may be lost. As a result, solutions to the difference 
equations may not be unique and may result in nonphysical spatial 
oscillations. Several different methods have been proposed to 
modify the CDS to guarantee positivity of coefficients. Probably the 
most popular of these schemes are the upwinding [50], hybrid [51], 
and power law [52]. These methods are well-behaved, but they 
cause high numerical diffusion. Recently, Runchal [53] proposed 
the CONDIF method which is designed to yield positive coefficients 
and, at the same time, retain the second-order accuracy of CDS. 
Although it was developed specifically to treat the cell-Reynolds- 
number problem in incompressible flows, there have recently been 
conjectures regarding its applicability to compressible flows with 
shocks, particularly its possible total variation diminishing (TVD) 
properties. In the sample problem examined below, it will be 
shown that this method is stable but that it is not TVD. Thus, it is 
not generally suitable for shock capturing. 



In this model problem, a nonlinear filter, based on a concept 
introduced by Engquist et al. [49], is used to study the dynamics of 
compressible vapor flow. The cell-Reynolds-number problem is 
involved as are shock capturing and shock reflections for 
intermediate and high Reynolds numbers. Other schemes 

(upwinding, power law, and CONDIF) will also be used to solve the 
problem, showing the numerical diffusion resulting from the use of 
these methods. 

4.2 SAMPLE PROBLEM 

The dynamics of compressible vapor flow are studied as a 
model problem. The vapor is assumed to be initially stagnant. A 
sudden increase in the input mass flux is imposed at a higher 
temperature. An example of this problem is the vapor flow 
dynamics of heat pipes. Fig. 4.1, in the transient mode. The 
governing equations are the usual continuity, momentum, and 
energy equations of compressible fluid flow. An equation of state 
relates the pressure to the temperature and density. These 
equations in normalized, one-dimensional (1-D) form are 
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VAPOR 





x = L 
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Fig. 4.1. One-dimensional vapor core model of a heat pipe 




dpT | dpuT d ( 1 
dt dx dx V RePr dx ) 


(4.3) 
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P = pT , 


(4.4) 
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and are solved on the domain (x,t) e (0,1) x (0 ,fy). In these 
equations, p is density, u is velocity, T is temperature, and p is 
pressure in nondimensional form. The nondimensional parameters 
are Reynolds number, Re, and Prandtl number, Pr, defined as 


Re = 


puD 



(4.5) 


where p is dynamic viscosity, k is thermal conductivity, C p is 
specific heat, and D is the pipe diameter. 


54 


I 



The initial and boundary conditions are 


p(x, 0) = 1, 


pu(x, 0) = 0, 


(4.6) 


T(x, 0) = T0(x ), 


and 

p(0, t) = p e , 


pu( 0, t ) = 1, 


(4.7) 


mo = T e , 


where p e and T e are, respectively, the density and temperature of 
the vapor at the evaporator of the heat pipe. 

4.3 NUMERICAL SOLUTION 

The governing Eqs. (4. 1 )-(4.3) are solved for p, pu, and T, 
respectively, with p recovered from the equation of state [Eq. 
(4.4)]. For the sake of analysis, the nonlinear partial-differential 
equations are cast in the general form 


(4.8) 
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d_ 

dt dx 


{^.rim.5. 

dx 


This equation is discretized on the staggered grid shown in Fig. 4.2. 
The temperature, density, and pressure are calculated at points 
indexed by i, and the mass flux, pu, is calculated at points indexed 
by /, where / = i + 1/2. 

Using the backward Euler algorithm for time differencing 
and centered differences in space, Eq. (4.8) is approximated by 


* + a 2,i?i - a 3,i 0 i + 1 = b i > 

where 

a l,i = D /-l +2 Pe t-^ ' 

a 3,i — D j (1 -J Pe { ), 


a 2.i ~ a l,i + a 3.i + ( u t~ + 


Ax 

At 


, „ Ax n 

b ; = Ax Si + <P- 

1 At 1 


with 




t-i 


Pi-i Ax 


= 


Pm** 


u, 

Pe ^D t ’ 


(4.9) 


(4.10) 


(4.11) 
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Fig. 4.2. The staggered grid and nodes 

I 

For stability and a unique solution, the coefficients of Eq. (4.9) 
should be positive. If the grid Peclet number, Pe is greater than 

I two in magnitude, some of the coefficients defined by Eq. (4.10) 

are not positive. 

I 4.3.1 Standard Schemes 

The coefficients of Eq. (4.9) are written in a general form as 

, a i,i = D /-i MPe^j) + max^M^.y) , 


a 3 i = D / A(Pe / ) + max(0, -u / ) , 

l 

Ax 

a 2,i = a l,i + a 3,i + ( u f u /-0 + ~J~ t * (4.12) 

I 

Expressions for A(Pe) for different schemes are listed in Table 4.1. 
Derivation of A(Pe) for upwinding and power-law schemes can be 


I 
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TABLE 4.1 


The Function A{Pe) for Different Schemes 


Scheme 

A (Pe) 

CDS 

1 - 0.5 1 Pe 1 

Upwinding 

1 


5 

Power law 

max[0,(l - 0.1 \Pe 1) ] 

CONDIF 

A (Pe t j) = 1 - max(0,^^ j/2) +R max (0,Pe t /2) 
A (Pe e ) = 1 - max(0, - Pe ^ 2) + -^-max(0, ]/2) 


found in Patankar [42] and for the CONDIF scheme in Runchal [53]. 
The parameter R appearing in Table I for the CONDIF scheme is 
defined as 


R = 


* 1+1 - 


(4.13) 


If there is a sharp variation in <J> then R may be quite large and, 
conversely, if there is only a small variation in <J>, R will be small. 
Therefore, a limit is imposed on R so that 







l/Rmax — ^ — Rmax 


(4.14) 


It can be seen, using the appropriate formulas from Table I, that 
Rmax controls the amount of numerical diffusion in the scheme. In 
particular, large values of R result in significant numerical 
diffusion while small values cause little diffusion. Typical values 
are in the range 4 < R ma x ^ 10. 

4.3.2 Filtering 

Engquist, et al. [49] introduced nonlinear filters for shock 
computation. These filters are used as post-processors in 
conjunction with standard finite-difference schemes. The 
difference equations are solved at every time level; then the 
solutions are processed with the filter to eliminate oscillations. 
The filters control the total variation (i.e., they possess the TVD 
property), produce sharp shocks (they are second order away from 
shocks), and they are simple to use. 

For the calculations, we use CDS and a nonlinear filter, as 
described below (for details, see Ref. [49]). The results are 
compared with those from the widely used upwinding, power-law, 
and CONDIF schemes. The filter produces more accurate results 
with no significant increase in the computation time for problems 


either experiencing cell-Reynolds-number difficulties or having 
shock-like solutions. 

The principles of the filtering technique are illustrated below 
using Fig. 4.3. Figure 4.3(a) shows the solution & n at time t n . 
Using CDS, the typical solution at time t n+1 , <p n+1 , is shown in Fig. 
4.3(b). This solution exhibits an overshoot, as would be expected. 
Results of filtering are displayed in Fig. 4.3(c) where it can be seen 
that the overshoot has been completely removed. The principles of 
the filtering process are the following: 


n+1 


(i) cp io is detected as an overshoot via the test 


,, , n+1 . , . n+1 .. . 

[(A.(p io )(A+(p in )]<0, 


r io 


(ii ) is decreased, 

n+1 

(iii) The new value of (p should not be less than the value 


- n+1 

of <Pi 0 -l * 


n+1 


n + 1 


(iv) (pi o + j is increased by the same amount as (p • is 


decreased in order to maintain conservation 
properties, and thus to obtain the right shock speed. 


The symbols A. and A + denote the backward and forward 
differences, respectively, and are defined as A± C> i = ±(O i:i:1 - <!>,). 
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Fig. 4.3. Filtering process applied to the CDS result 


The algorithm for this filtering and also for more complicated 
ones, can be found in Ref. [49], in which the nonlinearity and TVD 
properties of the filtering technique have been shown. 

4.4 RESULTS AND DISCUSSION 

The sample problem described above has been solved 
numerically using different schemes (centered differencing, first- 
order-upwinding, power-law, and CONDIF in conjunction with the 
nonlinear filter). Calculations have been carried out for high 

Reynolds number (Re = 10,000) to study the shock capturing 
ability of different schemes and also for low Reynolds number 
(Re - 20) to examine the cell-Reynolds-number problem. The 
calculations were performed on the Cray-2 supercomputer at the 
Lawrence Livermore National Laboratory. 
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In heat-pipe vapor-dynamics applications, the shock wave 
created in the evaporator travels towards the condensation region, 
where an outflow (smaller than the inflow at the evaporator) 
leaves the system, and the shock reflects back toward the 
evaporator. In this section, the reflection of a shock wave on a wall 
with zero permeability and with partial permeability has also been 
studied. 

4.4.1 Shock Capturing 

Figures 4.4 through 4.6 present the results for Re = 10,000. 
At this high Reynolds number, the diffusive terms are negligibly 
small, and the flow propagates mainly by convection. Figure 4.4 
compares the shock capturing ability of different schemes at two 
different times (f = 0.2 and 0.75). These calculations were 
performed with the number of grid points, N =41. Here, it is 
shown that the nonlinear filter leads to a sharp shock front, 
whereas the other schemes cause an overshoot and smearing of 
the front. The power-law-scheme results are basically the same 
results as those of upwinding for this case. This is due to the large- 
grid Peclet number. For the normalized velocity, length, and time, 
the shock front at t = 0.2 and t = 0.75 should be at x = 0.2 and 0.75, 
respectively. Although at early times all schemes calculate the 
shock at the right position, at later times ( t = 0.75) only filtering 
captures the shock front at the correct location. Comparing 
CONDIF and upwinding, it seems that the CONDIF scheme produces 
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Fig. 4.4 Mass flux profile at different times for different schemes 

a sharper front, but it causes more overshoot and does not catch 
the shock at the right position. 

In order to show how accurately the filtering captures the 
shock, calculations have been carried out for coarse and fine 
meshes. Figure 4.5 shows these results for the number of grid 
points N = 21, 41, and 81 at early ( t = 0.2) and later (f = 0.75) 
times. It is shown that the shock is captured within only two grid 
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Fig. 4.5 Mass flux profile at different times for different 

number of nodes 

spacings regardless of the number of grid points. Therefore, the 
shock-front width decreases significantly with a finer mesh. 

The total variation of the solution function, 0 , should not 
increase with time for a shock-capturing scheme with no 
oscillations or overshoot. The total variation at any time is defined 
by 
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N-h 

TV(0 n )= I 

i = / 


Then, the TVD property is expressed as 


TV(O n+1 ) < TV(0 n ) 


(4.15) 


(4.16) 


This criterion has been checked for the different schemes used in 
this study. Figure 4.6 shows that filtering possesses the TVD 
properties, whereas the upwinding and CONDIF methods do not 
exhibit TVD behavior. 

4 . 4.2 Cell-Reynolds-N umber Problem 

Nonlinear filtering has also been compared with several 
other schemes for treating the cell-ReynoIds-number problem. In 
this case the Reynolds number is relatively low, and the diffusion 
term is no longer negligible. However, for a coarse mesh, the grid 
Reynolds number may be greater than two, which causes 
oscillations. For Re = 20, the sample problem was solved for the 
mass flux, pu , on a fine mesh ( N = 81) using CDS. The CDS solution 
is compared with those of other schemes on a coarse mesh, N = 9, 
in Fig. 4.7 where it is seen that the CDS solution has an overshoot. 
Normally, when the cell-Reynolds-number is a problem, the CDS 
causes oscillations. However, as a result of the high compressibility 
of the flow in the sample problem, the pressure term in the 
momentum equation is of the same order as the other terms, 
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Time,t 


Fig. 4.6 Total variations of the mass flux for different schemes 

causing an overshoot due to forcing rather than oscillations. 
Although all of the schemes tested eliminate the overshoot, the 
filtering scheme produces the most accurate solution. To show the 
discrepancy from the "exact solution," the CDS solution on the fine 
mesh is taken as the reference solution, & r , and the L 2 -norm of 


the error has been calculated as 
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Fig. 4.7 Mass flux profiles for low Reynolds number 
for different schemes 



Table 4.II shows the £2 -errors for pu calculated at t = 0.5 for the 
various schemes considered. As can be seen, use of the nonlinear 
filter results in a significantly smaller error than any of the 
other methods. 
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TABLE 4. II 

L ? -Error for Different Schemes 



L 2 -norm 

CDS 

8.98 x 10' 2 

Upwinding 

6.37 x 10' 2 

Power law 

5.14 x 10' 2 

CONDIF 

3.06 x 10‘ 2 

Filtering 

1.49 x 10' 2 


4.4.3 Shock Reflection 

In this section we examine the use of nonlinear filtering in 
simulating reflected shocks on a partially blocked wall at x = L. 
When a shock reflects off a wall, the pressure on the wall increases 
to twice the value behind the shock [54]. This was verified in the 
computations by solving the equations for high Reynolds number, 
Re = 10 4 , with the boundary condition pu ( x =L ) = 0. In order to 
explore the shock reflection on a partially blocked end, the 
governing equations were first solved for low Reynolds number. 
Re = 10. In this case, the number of grid points was N - 21, and 
no filtering was needed. The physical process is initiated by a 
sudden inflow at x = 0. The mass flux profile has the same shape as 





I 


that shown in Fig. 4.7, and the corresponding wave front 
I propagates along the x -direction with time. When it reaches x = L, 

where only a fraction of the flow is extracted, the vapor is 
compressed and the pressure increases. This will cause a positive 
| pressure gradient along x, and the front reverses direction. The 

mass flux profiles along x at different times are shown in Fig. 4,8. 
When the flow front reaches x = 0, the vapor is compressed all 
I along x , and the inflow at x = 0 will cause the flow propagation 

towards the end x = L. When the flow front reaches the pipe's end 
for the second time, the increase in pressure causes more outflow, 
I and flow reflection takes place. The process of refill and reflection 

of mass flux continues, and the outflow increases with each 
reflection until the outflow is as high as the inflow. At this time 
I the process reaches a steady-state condition (Fig. 4.8.) 

The same vapor-flow transient behavior is found at higher 
Reynolds number. Figure 4.9 depicts the mass flow profiles at 
I different times for Re = 1000. In this case, the number of grid 

points is N = 41, and the calculated mass flow is filtered in two 

steps (i.e., forward and backward) at each time step. In the 
I forward step, the filtering process (Section 4.3.2) is used from 

x = 0 to the x -location of the shock front to eliminate the 
nonphysical oscillations. In the backward step, the mass flow 
I profile is filtered from x = L backward to the shock front. 

Referring to Figs. 4.8 and 4.9, the transient process reaches a 
steady-state condition after approximately five sets of refill and 


l 




I 


70 


I 




» 



I 


I 


Fig. 4.9 Shock reflections at high Reynolds number, Re - 1000 
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reflection of flow front, nearly independent of the input mass flux 
(or Reynolds number). 

4.5 CONCLUDING REMARKS 

The use of a nonlinear filter has been examined in the 
context of the dynamics of compressible vapor flow. The filter was 
used with both low and high Reynolds numbers to examine the 
cell-Reynolds-number problem as well as the shock-capturing 
ability of the technique. The resulting solutions were compared 
with those obtained by using standard schemes like upwinding, 
power-law, and CONDIF. The filtering was found to be quite simple 
to use, and is significantly more efficient than earlier methods with 
comparable ability. 

For low Reynolds numbers and a coarse grid (cell-Reynolds- 
number problem), the centered difference scheme produces an 
overshoot. Although all of the standard schemes mentioned above 
eliminate the overshoot, filtering produces the most accurate 
solution. Among these methods, upwinding produces the highest 
numerical diffusion, and the CONDIF method produces a more 
accurate solution than the upwinding and the power-law schemes. 

For high Reynolds numbers, the filtering technique is able to 
capture the shock in only two or three grid spacings regardless of 
the number of the grid points, whereas the other schemes cause an 
overshoot and smearing of the shock front. Comparing CONDIF and 


» 

upwinding, the former produces a sharper shock front, but it 
| causes more overshoot and the predicted wave speed is less 

accurate. 

Nonlinear filtering was also found suitable for calculation of 
| shock reflections. This aspect of the filtering procedure has not 

been tested in earlier work. When the wave front reaches a wall 
with partial permeability, the vapor compresses and the pressure 
I on the wall increases. This increase in pressure causes a backflow 

until the vapor is compressed over the entire length. 
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Chapter V 

STARTUP TRANSIENT OPERATION OF HEAT PIPES 


5.1 INTRODUCTION 

The transient behavior of vapor flow during the startup 
phase of heat pipe operation has been analyzed in order to 
examine the flow patterns for low- and high-input heat fluxes. 
Startup is a transient process through which the heat pipe starts 
its operation from a static condition and arrives at steady-state 
operation. The response of the vapor flow is an important design 
issue in dynamic thermal systems where a heat pipe is used to 
transfer heat following sudden increases in the heat load. 

During a startup transient, the working fluid of a heat pipe 
may initially be very cold or frozen in the wick. Here the liquid is 
assumed to be very cold and the vapor, which is 

thermodynamically in equilibrium with the liquid, is at very low 
pressure. The vapor core response to a sudden-input heat flux or a 
sudden change in the condenser temperature is the primary focus 
of this work. 

The vapor flow in heat pipes is a complicated problem 
because of the highly nonlinear nature of the governing equations 
and because of the inflow and outflow boundaries in the 



evaporator and condenser. Different approaches have been used to 

simplify the problem. In most previous work, the vapor flow is 

analyzed under steady-state conditions as a one-dimensional (ID) 
flow (e.g., Edwards and Marcus[12]) and a 2D flow [17-19, 24]. In 

studies of the dynamic behavior of heat pipes, the vapor flow is 

modeled either as a quasi-steady flow [38] or as a ID transient 
vapor flow along with friction coefficients approximated from a 
steady-state 2D model [37]. It has been shown by steady-state 2D 
analysis that the ID vapor-flow model is not able to accurately 
predict the axial heat and mass transfer and pressure drop [14]. 
Furthermore, 2D steady-state studies indicate that flow reversal 
takes place in the condensation section under high heat flux [41], 
It is therefore important to establish the conditions for this mode 
during transient operations. 

In a theoretical study using a ID transient model, Bystrov 
and Goncharov [31] showed that the time steps in the numerical 
calculations which are necessary to satisfy stability requirements 
are extremely small (order of 10' 7 ). The vapor dynamics of heat 
pipes was analyzed in Chapter Three, using a 2D transient model 
with the SIMPLER method. The primary advantage of SIMPLER is 
the staggered grid which makes the numerical scheme very stable 
and, consequently, the computational time very low. Although this 
method yielded interesting results, the resulting computational 
code is limited to low compressibility flows. However, in the 
startup mode of a heat pipe, there are regions of high as well as 
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low compressibility. Existing numerical codes for highly 
compressible flows, on the other hand, are only able to handle 
steady-state conditions because of their extremely small allowable 
time steps and, furthermore, they are not stable for low Mach- 
number problems. 

The objective of this chapter is to thoroughly investigate the 
dynamic behavior of the vapor flow in heat pipes during the 
startup transient phase of its operation. The fully compressible 
governing equations of vapor flow are numerically solved using a 
transient 2D model, discussed in Chapter Two. The numerical 
scheme is fully implicit using a staggered grid, and it is stable for 
all Mach numbers between zero and one, which is a suitable range 
for the startup transient mode. 

5.2 SOLUTION METHOD 

The complete set of governing equations of the vapor flow in 
heat pipes were derived in Chapter Two. These equations in their 
nondimensional form are Eqs. (2.11) through (2.15). The boundary 
conditions associated with these equations are given by Eqs. (2.17) 
to (2.22). The vapor is initially stagnant and at thermodynamic 
equilibrium with the liquid in the wick structure. When the initial 
pressure and temperature are very low, and a rather high input 
heat flux is applied in the evaporator, a shock wave is created in 
the evaporation zone. This is shown in Fig. 5.1 for a liquid sodium 
heat pipe. The vapor is at, initially, Po = 10"2 atm and T 0 = 800 K, 
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• Fig. 5.1. Vapor flow in a heat pipe at the initial stage of the startup 

and an input heat flux of 10 5 W/m 2 is applied at one end of the 

* pipe. The high gradients at the shock wave makes the SIMPLER 

method, used earlier, inappropriate for the startup transient 
analysis. In Chapter Four the one dimensional forms of the 

I governing equations were solved using standard methods such as 

upwinding and power law and a more recent scheme, CONDIF. It 
was shown that these schemes are well-behaved but that they 
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cause oscillations and overshoot at high compressibility and, 
except for CONDIF, they cause high numerical diffusion at low 
compressibility as well. The CONDIF scheme was also used by 
Issacci et al. [55] in a study of startup transient process in heat 
pipes. It was shown that for input heat fluxes more than 1 kW/m 2 
in an Na-filled heat pipe, a shock is created in the evaporation 
region. Therefore, the CONDIF scheme was found to be unsuitable 
for shock capturing. 

In Chapter Four it was shown that a centered difference 
scheme (CDS) used with nonlinear filtering yields a second order, 
stable solution in cell Reynolds-number problems and is able to 
capture a shock without oscillations. In this chapter, the same 
technique is used to solve the 2D form of the governing equations 
for a compressible vapor flow. 

The governing equations are nonlinear and coupled. An 
iterative method is used to solve each equation separately, i.e., 
equations (2.11) to (2.15) for p, pu, pv, T, and p, respectively. The 
SOR method is used to solve the differenced equations. Iterations 
on the equations were stopped after convergence of all the 
variables. The convergency was assumed to be reached when the 
relative change of all variables in one time step was less than 0.5%. 
This iteration method required the least storage and fewest 
calculations. 



5.3 RESULTS 


The boundary conditions of inflow and outflow at different 
locations on the same boundary are highly nonlinear. 
Furthermore, the exit point of the evaporator and the entrance 
point to the condenser, on the bottom boundary, are 
mathematically singular points. The boundary conditions at these 
points have sharp gradients which perturb the numerical scheme. 
In order to avoid this problem, the inflow and outflow mass fluxes 
are multiplied by the factor tanh X(x-x Q ). In the evaporator and in 
the condenser, x Q = 0 and 2, respectively. Depending on X, this 
factor makes the gradients smoother at the singular points. 

For the results shown in this chapter, the working fluid is 
liquid sodium at, initially, P o - 10*2 atm and T 0 = 800 K. The 
geometry dimensions are b - 5 cm, Le = L a = L c = b, and X = 5. The 
results are shown in the following sections for low- and high-input 
heat fluxes and for the pressure drops along different sections of 
the heat pipe. 


5.3.1 Error Analysis 

The calculational grid for each computation was chosen by 
inspection of the calculated L 2 - norm. The L 2 -norm is defined as 
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L 2 - norm = 


(5.1) 


where N x and Ny are numbers of grid points in x- and y- 
directions, respectively. Here, <p is one of the dependent variables 
(p, pu, pv, T, or P) and <p r is the calculated 0 on the finest mesh of 
calculations. Figure 5.2 shows the results of the error analysis for a 
sample calculation. With an increase in grid points, the L 2 -norm 
error decreases exponentially, however, the computational time 
increases (Fig. 5.2a). In this sample analysis, for grid points greater 
than 21, the error does not decrease by much, whereas the 
corresponding computational time increases significantly. 
Therefore, for this case, the optimized grid point is 21. 

The same error analysis is also carried out for time 
increments, At, using equation (5.1). In this case, <p r is the result of 
calculations according to the smallest At. Figure 5.2(b) shows the 
L 2 • norm error for different time increments. As expected, the 
error decreases as At decreases. The sample error analysis, 
described above, was carried out for different input conditions in 
order to choose the corresponding optimized grid points and time 
increments. For all the calculations presented below, At = 10' 2 to 
10 - 3 . 

5.3.2 Low-Input Heat Flux 

Figure 5.3 shows the flow patterns in the vapor core for a 
low-input heat flux, q 0 =10 5 W/m 2 . The Reynolds number based 

on the vapor thickness is Re = 100 and the computations were 
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Fig. 5.2 L 2 -norm error and computational time 
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done for an optimized grid of 21 x 61. The transient development 
of the vertical mass flux, pv, at different times is shown in Fig. 5.3. 
At t= 0, the vapor is stagnant. Applying the input heat flux, 
evaporation takes place. Since the input heat flux is relatively 
high, a shock wave is created (Fig. 5.3a). The vapor flow develops 
above the evaporator and in the adiabatic region and the shock 
travels above the evaporator until it hits the upper boundary (Fig. 
5.3b) where the vertical mass flux is blocked, pv = 0. At this point, 
the vapor is compressed and the vapor pressure increases which 
causes a positive vertical pressure gradient and the shock reflects 
back (Fig. 5.3c). The figure also shows reverse flow in the top left 
comer which is caused by shock reflection. When the reflected 
shock reaches the lower boundary (Fig. 5.3d) and the vapor is 
compressed within the evaporation region, the vertical pressure 
gradient is reduced and the inflow mass flux tries again to fill the 
evaporation region (Fig. 5.3e). As the vapor fills the evaporation 
region (Fig. 5.3f), another shock reflection occurs and the flow 
reverses in the top left corner (Fig. 5.3g). Comparing Figs. 5.3d and 
5.3g, the reverse flow after the second reflection is weaker than 
that in the first reflection. About eight to ten cycles of refill and 
reflection of mass flow are observed in the evaporation region 
before the process becomes steady (Fig. 5.3i). At each cycle the 
reverse flow is weaker and after four cycles there is no reverse 
flow. The vertical gradient in the mass flux, pv, also decreases with 
each cycle. 
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Fig. 5.3 (Con’d). Transient devi 
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The transient flow pattern in the condensation region, 
2<x/b ^3, is also illustrated in Fig. 5.3. At early stages of the 
startup transient process, the vapor is at the same temperature as 
the liquid and, therefore, there is no condensation (Figs. 5.3a and 
5.3b). When vapor at the higher temperature flows from the 
evaporation to the condensation region, the vapor temperature in 
the condenser increases. A temperature gradient is then 
established along the liquid layer since the temperature at the 
outer surface of the liquid layer is kept fixed, and condensation 
takes place (Fig. 5.3c). The vapor temperature in the condenser 
increases with time and, consequently, the condensation rate 
increases. At steady state, shown in Fig. 5.3i, the mass flow rate 
through the evaporator equals the mass flow rate out of the 
condenser. 

5.3.3 High-Input Heat Flux 

In this section the transient flow patterns are shown for a 
relatively high heat flux, q Q = 10 6 W/m 2 . The corresponding 

Reynolds number in this case is Re = 1000 and the optimized grid 
is 41 x 121. Figure 5.4 shows the transient development of the 
vertical flow, pv, at different times. In the evaporator region, 
multiple shock reflection and refill of vapor is observed to be the 
same as that for the low-input heat flux discussed above. In this 
case, however, the reversed flow caused by the shock reflection 
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Fig. 5.4 (Con'd). Transient development of the vertical mass flux. 

High-input heat flux 
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develops with time and is sustained as the process reaches steady- 
state operation. 

At high-input heat flux the flow pattern in the adiabatic 
region, 1 <x/b £ 2, is significantly different from that for lower 
heat fluxes. Figure 5.4d shows that when the shock reflection 
reaches the lower boundary of the evaporator, circulation is 
initiated in the adiabatic region. This circulation ceases when the 
vapor flow refills the evaporator region (Fig. 5.4e). However, the 
second reflection will cause a stronger circulation in the adiabatic 
region, as shown in Fig. 5.4f. Figures 5.4g through 5.4/ 
demonstrate the subsequent sets of reflection and refill of the 
vapor flow with circulation strength increasing with time. When 
the process reaches steady-state conditions, circulation occupies a 
significant portion of the region. Circulation in the adiabatic region 
was also observed in Chapter Three in the analysis of the 
operating transient mode of heat pipes. 

In order to show the overall picture of flow pattern 
development. Fig. 5.5 was prepared to depict the vapor flow fields 
at different times. In this figure the arrow indicates the flow 
direction and the tail indicates the magnitude of the flow. Figure 
5.5a shows the flow field at an early stage of the transient process 
after the first shock reflection. Here the flow reversal is shown to 
be in the left top corner of the evaporation region. The shock 
reflection initiates flow circulation at the entrance to the adiabatic 
region shown in Fig. 5.5b. Vortex formation in the adiabatic region 
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Fig. 5.5. Vapor flow fields at different time levels 
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Fig. 5.5 (Con'd). Vapor flow fields at different time levels 
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i develops with time to a stronger flow circulation (Fig. 5.5c) and 

another vortex can be seen in the left corner of the evaporation 
region. In the flow field at steady state (Fig. 5.5d), it is seen that 

0 flow circulation in the adiabatic region develops to two vortices. 

One of the limitations on the heat pipe operation is the sonic 
limit. When the vapor velocity reaches the speed of sound the 

1 vapor flow is choked, which limits the heat flux through the pipe. 
Therefore, the maximum Mach number in the vapor flow should 
be equal or less than one. Figure 5.6 shows the development of the 

I Mach number profile with time along the pipe, for high-input heat 

flux ( q Q = 10 6 W/m 2 , corresponding to the results presented in this 

section). It is shown that the Mach number is below one at the 
I initial stage of the transient process. The vapor flow develops in 

the evaporation region and the Mach number increases with time. 
The vapor velocity reaches the speed of sound at the exit of the 
I evaporator, where the Mach number is equal to one. 

5.3.4 Pressure and Pressure Drop 

I One of the major issues in designing a heat pipe is to 

correctly calculate the pressure drop along the vapor core in 
different sections of the pipe. The correlations widely used in 
I these calculations are based on simple ID steady-state analysis. In 

this section, the pressure drops in different sections of the heat 
pipe are shown as functions of time. The calculated pressure drops 
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Fig. 5.6. Mach number profiles at different times; 
high-input heat flux 

for high- and low-input (q 0 = 10 6 W/m 2 , and 10 4 W/m^, 

respectively) heat fluxes are shown in Fig. 5.7 for different times 
during the transient phase. Solid lines show the maximum 
pressure drop in a specific region of the pipe. Averaging the 
pressure along the vertical cross section of the pipe and then 
calculating the maximum-averaged pressure drop yields 
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Fig. 5.7 Calculated pressure drops for (a) high- and (b) low-input 

heat fluxes ( and ••• = maximum and maximum- 

averaged pressure drops, respectively) 
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interesting results which are shown (dotted lines) for different 
times in different sections. Both the maximum and the maximum- 
averaged pressure drops display nearly periodic oscillation until 
they converge to a constant value at steady state. Pressure 

oscillations are caused by multiple shock-wave reflections in the 
evaporation region (Figs. 5.3 and 5.4) and at each shock reflection, 
the local pressure increases significantly. Consequently, a large 
pressure difference is created along the pipe. 

Another interesting result which can be seen on Fig. 5.7 is 
that the pressure drop is significantly higher during the initial 
period of the startup transient of heat pipe operation. The mean 
value of the pressure drop decreases with time, and at steady 

state its value is less than 10% of the initial pressure drop. 

Therefore, if a heat pipe is designed to operate in the transient 
mode, the use of steady-state correlations may result in significant 
errors. 

Tien and Rohani [18] showed that the cross-sectional 

averaged pressure along a heat pipe is greater than the pressure 
calculated by Bankston and Smith [15] using a ID analysis. Here 
we compare the cross-sectional averaged pressure with the 
calculated 2D pressure by comparing the maximum pressure drops 
in both cases. Although the maximum-averaged pressure drops 
(dotted lines on Fig. 5.7) converge to the maximum pressure drop 
as the process approaches steady state, it is significantly lower 
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during the initial transient time. This emphasizes the need for a 2D 
vapor flow model in an overall heat pipe transient analysis. 

The pressure drop along the pipe, as mentioned before, plays 
an important role in the design of a heat pipe. In addition, the 
pressure in the evaporator may cause limitations on heat pipe 
operations. When a shock wave reflects in the evaporation region 
and the local pressure increases, the increase in the overall 
pressure of the evaporator may affect the evaporation rate. In 
order to study this aspect, the averaged (solid lines) and the 
saturated (dotted lines) vapor pressures in the liquid-vapor 
interface in the evaporator are shown as functions of time in Fig. 
5.8. Here again, the periodic behavior of pressure is evident. 
During early transient times, the averaged pressure in the 
evaporator is sometimes larger than the saturated pressure. This 
will affect the evaporation process by decreasing the evaporation 
rate. Furthermore, whenever, the vapor pressure is greater than 
the saturated pressure, condensation takes place by nucleation of 
mists in the vapor. In this case the flow is not homogeneous 
anymore and a two-phase flow analysis should be incorporated. 

5.3.5 Friction Factor 

The total pressure drop in the vapor core is needed in 
designing a heat pipe. The averaged friction factor along the pipe is 
calculated by 
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Fig. 5.8 Averaged ( ) and saturated (•■•) vapor pressures in the 

evaporator 
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Figure 5.9 shows the ratio of the averaged friction factor and its 
value at steady state, C f (t)/C f (°°), for low- (Fig. 5.9a) and high- (Fig. 




C f (t)/C,(~) C,(t)/c,( 




Fig. 5.9. Friction factors, a) Low-input heat flux, 
b) high-input heat flux 




5.9b) input heat fluxes. The friction factor oscillates periodically 
because of the shock reflections, and converges to a constant value 
at steady state. Figure 5.9a shows that for low-input heat flux (q" = 
10 4 W/m 2 ), the maximum-friction factor during the transient 
phase, is about %50 more than its value at steady state. For high- 
input heat flux (q» = 10 6 W/m 2 ), the friction factor is about 35 
times more than its steady-state value (Fig. 5.9b). This implies that 
in designing a heat pipe for transient operations, the steady-state 
correlations may cause a significant error. 

5.4 CONCLUSIONS 

A nonlinear filtering technique has been used to analyze the 
startup vapor dynamics in heat pipes. Large variations in vapor 
compressibility during the transient phase render it inappropriate 
to use the SIMPLER scheme, as well as standard schemes for the 
disretization of the advective terms. 

For a high heat flux, the startup transient phase involves 
multiple shock reflections from the line of symmetry in the 
evaporator region. Each shock reflection causes a significant 
increase in the local pressure and a large pressure drop along the 
heat pipe. Furthermore, shock reflections cause flow reversal in 
the evaporation region and flow circulations in the adiabatic 
region. The vapor vortex formation in the evaporator is only 
transient and is observed to disappear as the dynamics approach 
steady state. However the circulation in the adiabatic region 


grows with time, and in the steady-state condition the circulation 
occupies a significant portion of the region. 

The maximum- and maximum-averaged pressure drops in 
different sections of the heat pipe oscillate periodically with time 
because of the multiple shock reflections in the evaporator. The 
pressure drops converge to a constant value at steady state. 
However, they are significantly higher than their steady-state 
value at the initiation of the startup transient. 

The pressure increase due to shock reflection may cause the 
overall pressure of the evaporator to be greater than the saturated 
vapor pressure. When this occurs, condensation takes place by 
mist nucleation in the vapor. Therefore, the condensation rate is 
expected to decrease. 

The time for the vapor core to reach a steady-state condition 
depends on the input heat flux, the heat pipe geometry, the 
working fluid, and the condenser conditions. However, the vapor 
transient time is on the order of seconds. Depending on the time 
constant for the overall system, the vapor transient time may be 
very short. Therefore, the vapor core may be assumed quasi- 
steady in transient analysis of a heat pipe operation. 



Chapter VI 

EXPERIMENTS: NEEDS AND DIFFICULTIES 


6.1 INTRODUCTION 

Transient behavior of the heat pipe internal flows needs to 
be examined experimentally to understand and properly model 
the physical phenomena governing their behavior. The 
experimental data are needed to be compared with and verify the 
theoretical results. Experimental studies of heat pipe dynamics are 
very limited in the literature. Most of the experiments were 
carried out to measure only the maximum possible input heat flux 
and/or the temperature of the heat pipe external surface [27, 28]. 
These data are limited to few working fluids and operational 
conditions. Experimental data for the vapor flow in heat pipes are 
not available because of the complexity and sensitivity of the flow 
to external objects. The standard measurement methods like hot 
wire anemometry or thermocouples cannot be used to measure the 
vapor velocity and temperature. These measurement tools will 
disturb the energy transport process in the pipe and the results 
will not be reliable for comparison with the analysis. 



The following is description of an experimental heat pipe, 
designed for transient vapor flow measurements. Then the optical 
measurement methods, which are non-invasive and unique for 
heat pipe internal flow measurements, are discussed. Experimental 
studies on the dryout and rewet of the heat pipe wick structure 
are also included. 

6.2 EXPERIMENTAL APPARATUS 

Figure 6.1 shows a schematic diagram of the experimental 
apparatus. A rectangular pipe has been designed and built to 
simulate a heat pipe. The bottom of the pipe is a plane wide and 
thin channel which contains the working fluid. Two layers of 
stainless steel mesh (300 series, 400 mesh) are fixed on top of the 
liquid channel to serve as the wick structure. The axial sides of the 
pipe are made of glass to permit flow visualization. The top of the 
pipe is insulated. One end of the bottom plate is heated by patch 
heaters which supply up to 1.55 W/cm 2 . The other end is cooled 
by air flow. The middle region of the bottom plate is insulated to 
approximate the adiabatic conditions. 

Freon 11 and freon 113 are chosen as working fluids. These 
fluids have low latent heats and low boiling points and are easy to 
handle. The low latent heat allows one to induce a liquid-vapor 
phase transition at low applied heat fluxes. A scaling analysis 
might be used to generalize the results to other fluids. 
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Fig. 6.1 Experimental rectangular heat pipe 

The axial and vertical velocity components are measured 
using Laser Doppler Anemometry. The temperature and 
concentration profiles in the vapor core are measured 
simultaneously using twin beam laser holographic interferometry. 
These methods are described briefly in the following. 

6.3 OPTICAL METHODS 

Optical methods have been used for many years in heat and 
mass transfer measurements due their unique advantages. These 
methods do not disturb the examined process and the information 
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of the whole field is recorded on photographs, which has the 
advantage on the methods of point by point measurements. 
Furthermore, the measurements are inertialess and therefore, very 
fast phenomena can be investigated [56]. 

In this section, the use of optical methods in heat pipe 
experiments are discussed. The laser Doppler anemometry is used 
for measuring the axial and vertical velocity components of vapor 
flow. The advantages of the twin beam holographic interferometry 
are discussed as a means of measuring both the temperature and 
concentration distributions in the vapor core, simultaneously. 

6.3.1 Holographic Interferometry 

Interferometric methods are based on differences in lengths 
of the optical paths. Among these methods, holographic 
interferometry is the cheapest and easiest to construct with a 
comparative accuracy [56]. However, a highly coherent light is 
needed. The single beam holographic interferometry is first 
illustrated. Then it will be shown how this method can be modified 
to twin beam interferometry. 

Sin gle Beam Interferometry 

Figure 6.2 is a schematic diagram of a single beam 
interferometer. The light source is a laser. The laser beam is split 
into an object beam and a reference beam. Each beam is passed 
through a microscope objective and a collimating lens and 
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Fig. 6.2 Schematic of single beam holographic interferometer 
I 1. Helium-Neon laser, 2. Argon laser, 3. variable beam 

splitter, 4. object beam, 5. reference beam, 6. micro- 
scope, 7. pinhole, 8. collimating lens, 9. test section, 

I 10. holographic plate, 11. lens, 12. camera, 13. air- 

suspension table, 14. shutter 

I expanded to a parallel beam. The reference beam falls directly on 

the hologram plate, whereas the object beam passes through the 
experimental apparatus and then intersects the reference beam on 
! the hologram plate. 

104 


l 




I 


In real-time holographic interferometry the photographic 
I plate is exposed before the experiment starts, and the comparison 

beam is recorded. The photograph plate is developed and fixed 

while it stays in its place. The comparison beam is then 
I reconstructed by illuminating the hologram with the reference 

beam. After the experiment starts a phase change occurs in the 
object beam and illumination of this beam on the hologram 

I produces interference fringes. Changes in the interference pattern 

due to changes in the process can be continuously observed or 
photographed. 

I Interferometric methods are based on the measurement of 

variations of the refractive index in the test section. These 
variations may be caused either by temperature, concentration or 
I pressure gradient fields. Then any of these fields alone can be 

determined by evaluation of the interferogram. However, if the 
refractive index is changed due to temperature and concentration 
l changes simultaneously, an interferogram produced by a single 

beam interferometer may not be adequate to evaluate the 
temperature and/or concentration fields. This difficulty can be 
I overcome by using two laser beams with two different 

wavelengths. This method is called twin beam (or two-wavelength) 
holographic interferometry and was first proposed by Mayinger 
I and Panknin [57]. 
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I The first accurate results were obtained by Panknin [58]. 

Using twin beam interferometry he found good results for 

measurement of the temperature and concentration in 
I simultaneous heat and mass transfer along a heated vertical plate 

and a horizontal cylinder with sublimation and the temperature 
and fuel distributions in a flame. 

I Figure 6.3 shows the required experimental arrangement. 

The twin beam interferometer is very similar to that of the single 
beam shown in Fig. 6.2. Here, two lasers, an Argon laser with 
I wavelength Ay =457.9 nm and a Helium-Neon laser with 

wavelength X k =632.8 nm, are used as light sources. The two 
beams lie on a horizontal plane and are perpendicular to each 
| other when they intersect on the beam splitter. Therefore, both 

object and reference beams consist of two different wavelengths. 

The exposure technique in twin beam interferometry is 
I similar to that of the single beam. Before the experiment is started, 

the comparison beam, which now consists of two wavelengths, is 
recorded on a photograph plate. The photograph plate is developed 
I and fixed. Now, if the laser beam with the wavelength, say, A j , is 

blocked, both reference and object beams will have wavelength A*. 
Then by illuminating the reference beam on the hologram, the 
I comparison beam is reconstructed. The object beam which passes 

through the test section is illuminated on the hologram and 
produces an interference pattern for wavelength X k . By shining the 
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Fig. 6.3 Schematic of twin beam holographic interferometer 
1. Helium-Neon laser, 2. Argon laser, 3. variable beam 
splitter, 4. object beam, 5. reference beam, 6. micro- 
scope, 7. pinhole, 8. collimating lens, 9. test section, 

10. holographic plate, 11. lens, 12. camera, 13. air- 
suspension table, 14. shutter 

laser beam of wavelength Xj , another interference pattern is 
produced for wavelength Xj. 
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The main advantages of this technique is that only one plate 
is developed which includes information from both wavelengths, 
and two interference patterns are produced by shining each laser 
beam, separately. Two sets of fringes are, then, evaluated to 
determine the temperature and concentration profiles for the 
tested medium. The evaluation method of fringes is discussed in 
the following section. 

6.3.2 Theory of Twin Beam Interferometry 

The following assumptions are made in development of the 
governing equations [38]: 

1. The optical system is perfect, the experimental setup is 
mechanically stable and the lasers are ideal. 

2. The object beams with the wavelengths Xj and X k are ideal 

parallel waves. 

3. The variation of the refractive index is only two dimensional 

and no variation in the beam direction. 

4. No reflection of beam due to gradients of the refractive index. 

5. Refractive index is constant outside of the test section and 
during recording of the comparison beams. 

6. The holographic construction is perfect. 

With the above assumptions the interference pattern 
depends only on the variation of the refractive index. In 
holographic interferometry the object beam, passing through the 
test section at different times, is superimposed on the comparison 
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beam, and therefore, reveal the difference in optical pathlengths in 
two exposures. Expressed in multiple S of the wavelength A, this 
difference is written as 


S(x,y,A)A = l[n(x,y,X)~ rtj 


( 6 . 1 ) 


where / is the width of the test section, in which the refractive 
index varies because of temperature and concentration gradients. 
The interference pattern shows the change of the refractive index 
between the recorded comparison beam (constant temperature T m 
and constant concentration C m , yield a constant refractive n_) and 
the "measurement" beam. 

An extinction of light occurs for 


1 ' 2 2 2 


(dark fringes). 


( 6 . 2 ) 


and an amplification for 

|S| = 1,2,3... (bright fringes). (6.3) 

The interference fringes are points of the same refractive index 
change. With temperature and concentration gradients the 
refractive index variation depends on both gradients. Therefore, a 
relationship has to be found between the refractive index change 
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due to the temperature gradient and that due to the concentration 
gradient. 

The molar refractivity N(X), is related to the refractive index 
by the Lorentz-Lorenz equation as follows 


N(X) = 


n(X ) 2 - 1 M_ 
n(X ) 2 + 2 X p 


(6.4) 


where M is the molecular weight and p is the density. Fortunately 
the molar refractivity does not depend on temperature and 

pressure, but it varies with wavelength. This is the basis for the 

I 

twin beam interferometry. 

For gases with rt~U Eq. (6.4) is approximated quite 

accurately by the Gladstone-Dale equation 

I 

N(X) = -[n(X)-l\. (6-5) 

3 P 

I 

For the ideal gases the density, p, can be replaced by the total 
pressure, p, and the temperature, T, to get 

1 

N (X) = —[n(X)~ 1], ( 6 - 6 > 

3 P 

S where R is the gas constant. The molar refractivity, N(X ), for a 

mixture of two components is given by 
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N(X) = C a N a (X) + C b N b (X) t 


(6.7) 
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with C a + C b = 1. In the above equation C a and C b are the 
concentration of the components a and b in the mixture and N a (X) 
and N b ( X) are the molar refractivity of the components in their 
pure form. 

Combining Eqs. (6.1), (6.6), and (6.7) yields 

S(x,y.K).X = ^{^^[/V.(A) + C > U, 5 ')(iV l (A)-JV.a))]-i-W a (A)| i 

( 6 . 8 ) 


with C a = l-C 6 , C am = 1 and C a- = 1 (only component a presents in 
recording the comparison wave). For two different wavelengths Xj 
and X k we get two equations from which the temperature and 
concentration distributions will be evaluated. These equations are 
the following 


S(x,y, X t ) 


A N(X,) 


— S(x,y,X k ) 


■k 


T(x,y)-T„3pl 


r(x,y)T. 2 R 


A N(X k ) 

Hg(Xj) N a (X k ) ' 
AN{Xj) A N(X k ) 


(6.9) 




C(x,y) CS 

3 pi 

'AN(Xj) AN(X k ) 

T(x,y ) 7. 

' 2 R 

_N a (Xj) N a (X k )_ 


1 1 1 
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where AN(X ) = N a (X ) - N b (X ) . 
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Figure 6.4 demonstrates how to evaluate the temperature 
profile, T(y), and the concentration profile, C(y), from two 
interferograms of a boundary layer produced by two different 
wavelengths, Xj and X k [58}. For a specified direction along the 
boundary layer, x= constant. Fig. 6.4b shows Sj, S k , Sj. Xj and S k .X k 
as functions of y, where Sj = S(x =const., y, Xj) and Sj = S(x =const., 
y, Xj). According to Eqs. (6.9) and (6.10) the temperature and 
concentration profiles at a location are found by the following 
equations 
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T(y)-T m ~ 


C(y)-C. 


S -X 
1 1 

s t x, 

ANj 

AN k ’ 

Sj-Xj 

S k - X k 

N < 

N k ' 


( 6 . 11 ) 


( 6 . 12 ) 


where AT = N a (Xj) and N k = N a (X k ). Figures 6.4c and 6.4e show the 
variations of S-X/AN and SX/N for Xj and X k . Then using these 
distributions, Eqs. (6.11) and (6.12) yield the temperature and 
concentration profiles for that location. These profiles are shown in 
Figs. 6.4d and 6.4f. 
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(a) 


(b) (c) 



(d) 


(e) (0 


I Fig. 6.4 Evaluation of interferograms [58] 

a) interferograms for Ay and X k 

b) phase lag S • A, 5 

I c), e) difference of modified phase lags 

d), f) temperature and concentration profiles 

I 6.3.3 Laser Doppler Anemometry 

The laser Doppler anemometry allows the measurement of 
the local, instantaneous velocity of tracer particles suspended in a 
1 flow without disturbing the flow. Thus, appropriate particles must 

exist in the fluid and the relationship between the particles and 
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fluid velocity must be known. The task of a laser Doppler 
anemometer falls into three parts. 


• Generating two coherent laser beams which interfere in a 

probe volume 

• Detecting the light intensity scattered by tracer particles 

crossing the measuring control volume 

• Analyzing the Doppler frequency and calculating the 

velocity of the particles 

Figure 6.5 shows the required components of a laser Doppler 
anemometer and the following is a brief description of the basic 
principles of laser Doppler anemometry which are carefully 
outlined in Durst et al. [59]. 

Two coherent laser beams with wavelength A and a diameter 
d form an ellipsoidal volume at the beam intersection point (see 
Figs. 6.6 and 6.7). The size of the measuring control volume is 
given by 


5 x = d, 



(6.13) 


where <t> is the angle between the two beams and 8 X ,8 y and 8 Z are 
the axes of the ellipsoid. Figure 6.7 shows the principles of the 
"interference fringe" model which was proposed by Rudd [60]. 

Two coherent laser beams having plane wave fronts intersect 
at an angle, <p. This yields a pattern of plane interference fringes as 
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Fig. 6.5 LDA optics from DANTEC manual 

1. beam splitter 2. bragg cell unit 3. beam displacer 
4. beam expander 5. front lens 6. pinhole unit 
7. photomultiplier optics 8. photomultiplier 
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shown in Fig. 6.7. The fringe spacing , Az, is proportional to the 
wavelength. A, of the light and inversely proportional to the half 
angle between two beams 


A z = 


A 

2sin{<f> / 2) 


(6.14) 


As a particle moves through the region of interference fringes it 
will scatter light whose intensity will vary according to the light 
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Fig. 6.6 Measuring control volume 
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intensity variation inside the bisector of the two beams. These 
light intensity variations have the Doppler frequency v D given by 

U 2Usin(<p/2 ) 

- — • (6 - 15> 


where U is the particle velocity component perpendicular to the 
fringes. 

The scattered light is detected by a photomultiplier and the 
velocity component normal to the interference fringes is 
determined by analyzing the Doppler frequency. This frequency is 
ambiguous because it does not indicate in which direction a 


I 
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Fig. 6.7 Interference fringes in the control volume 


I particle moves. This lack of information can yield large errors if 

the flow direction is not known (e.g. in case of recirculating flow) 
or if the mean velocity is near-zero and the fluctuations are larger 
I than the mean value in turbulent flow. 

Frequency shifting allows the 180 degree ambiguity in the 
direction of a measured velocity to be resolved. The frequency 
I shift is obtained when the two intersecting laser-beams have 

different frequencies. The presence of a frequency difference, v 5 ,, 
results in a movement of the fringe pattern with velocity U s 
I 

U s = -v s Az . (6.16) 
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A frequency shift, applied to one of the laser beams, will result in 
an increase or decrease in the measured velocity of a particle 
depending on the direction of the movement. In this case the 
velocity can be determined clearly from the detected Doppler 
frequency, v D ,. If the velocity of the fringes, caused by the shift 

frequency, is greater than the particle velocity, then 


U = ( V D~ V s) 


2sin(<p 12) 


( 6 . 17 ) 


A lower Doppler-frequency than the shift frequency indicates that 
a particle moves in the same direction as the fringe pattern. If a 
particle moves in the opposite direction, the Doppler frequency 
will be larger than the shift frequency. 

Laser Doppler anemometry depends on signals from particles 
suspended in the flow. Thus, it is important to know how well the 
motion of the particles represent the motion of the fluid. In our 
case there are two reasons for velocity differences: 

• The particle terminal velocity is different from the vapor 
velocity because of greater gravitational forces on the particle 
as compared to the vapor. 

• Drag forces caused by accelerations and velocity fluctuations 
can result in particle/vapor mismatch. 
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According to Tchen [61] the terminal velocity, C f , of a particle 
in a fluid with the constant velocity Uj- parallel to the gravity g is 




where 




(6.18) 


(6.19) 


The velocity difference C f between the particle velocity, U , and 
the fluid velocity, U ^ , depends on the density difference, {p p -p f ), 
the viscosity of the fluid, p.^ and the squared radius of the 
particles, r 2 . For oil droplets (r,»l pm) in gas, at standard state, 
C, => 1.5x1 O' 4 m/s. 

According to Dring [62] the velocity difference between a 
particle and fluid caused by high-frequency velocity oscillations, 
depends on the Stokes number 


18p f 


( 6 . 20 ) 


For Stokes number St < 0.14 the particles follow an oscillation with 
the frequency o) with less than 1% velocity difference. For oil 
droplets in a gas, at standard states, ©^ = 12 kHz. 
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6.4 PRELIMINARY RESULTS 

Experiments have been carried out to visualize the vapor 
flow patterns in evaporation, adiabatic and condensation regions of 
the heat pipe after a sudden increase in input heat flux has been 
applied. A boundary layer of vapor is established very quickly on 
top of the wick surface all over the heat pipe. The thickness of this 
layer decreases from the evaporation section to the condensation 
region. The vapor boundary layer grows quickly and fills the 
whole space. Depending on the input heat flux, the whole process 
takes place in a few seconds. For low input heat flux, the vapor 
flow is found to be smooth and laminar. For higher heat fluxes the 
vapor boundary layer grows faster. If the input heat was high 
enough a few regular vortices were observed on top of the 
evaporation part and they expanded to the entrance region of the 
adiabatic part. The vapor flow in the adiabatic and condensation 
regions under a high heat flux are still laminar but reverse flow on 
top of the condenser is suspected. 

6.5 DRYOUT AND REWET OF THE WICK STRUCTURE 

A major limitation in the operation of heat pipes is wick 
dryout. The dry region of the wick will heat up, maybe even 
overheat, and therefore be difficult to rewet. The forces resisting 
dryout and augmenting rewet need to be determined before one 
can model a heat pipe to study its operational characteristics. Our 
search for data led us to conclude that it would be instructive to 
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first examine dryout and rewet of a thin film of fluid flowing over 
a flat, inclined plate, heated from below. 

The stability of thin films flowing over both heated and 
unheated surfaces has been extensively studied. Yih [63] and 
Benjamin [64] analytically studied the purely fluid mechanical 
problem in the absence of any heat load. The neutral-stability 
curves obtained by Benjamin [64] show that free-surface flows 
down a vertical plane are unstable for all Reynolds numbers. 
Later, Bankoff [65], starting with a perturbation procedure 
provided by Yih [66], studied the combined fluid mechanical and 
heat transfer problem. In his study evaporation from the surface 
of a falling liquid film was found to have a destabilizing effect. 

The stability of a dry spot, once it occurs, has been analyzed 
by several investigators. Murgatroyd [67], Zuber and Staub [68], 
and McPherson and Murgatroyd [69] predicted conditions under 
which a dry spot becomes stable. This was usually achieved by 
means of a force balance in which the pressure developed at the 
upstream stagnation point of a dryout region, the thermocapillary 
force developed as a result of the variation of surface tension with 
surface temperature, and other forces were taken into account. 

When a liquid film becomes ultra thin, long range molecular 
forces become significant and must be integrated with continuum 
mechanics in an analysis of dryout or rewet. Bankoff [70] studied 
ultra thin films analytically and found that nonlinear wave theory. 


bifurcation theory, and asymptotics all contribute to our 
understanding of the dynamics of dryout and rewet. 

Experimental studies of dryout and rewet in thin films 
flowing over a flat plate are rare. Orell and Bankoff [71] studied 
the formation of a dry spot in a non-boiling, thin film of ethanol on 
a horizontal surface subject to a slowly increasing heat flux. It was 
found that noticeable thinning of the film occurred prior to rupture 
and that the heat flux at which rewet occurred was, on the 
average, 30% less than the heat flux needed for dryout at the same 
film thickness and flow rate. It was noted that the flow rate had a 
significant effect on the heat flux required for dryout, but this 
trend was not formally investigated. 

In this section we focus on the effects of film thickness, flow 
rate, and inclination angle on the minimum heat flux required for 
dryout and the maximum heat flux allowed for rewet. 

6.5.1 Experimental Setup 

Fig. 6.8 shows the schematics of the experimental setup. The 
apparatus consists mainly of a flat copper plate whose dimensions 
are 35.5 cm x 15.2 cm x 0.24 cm. By adjusting the stopcock on an 
Oppenheimer tube mounted above the upper reservoir, we can 
control the flow rate through the knife edge gap. This gap is set 
between 0.025 mm and 0.080 mm depending on the fluid viscosity 
and the desired film thickness. 



Fig. 6.8 Schematic of film dryout experiment 

The last 11.4 cm of the plate constitutes the heated zone. 
Heat flux is provided by four 12.7 cm x 2.5 cm patch heaters which 
are capable of delivering 1.55 W/cm 2 at 115V. The plate 
temperature in this region is monitored using 15 copper- 
constantin thermocouples mounted 0.5 mm from the upper surface 
of the plate. 

The thickness of the films is measured using a micrometer 
head with a 30° point. This micrometer is mounted just above the 
surface of the fluid in the center of the heated region. Film 
thickness can be measured to within ±0.02 mm in this fashion. 
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6.5.2 Preliminary Results 

The fluids used to form thin films exhibited very different 
physical characteristics. During the course of experimentation, 
three general mechanisms of dryout were observed. The first of 
these is dryout due to rupture. This mechanism occurred when 
water was used as the working fluid. No thinning of the film is 
observed when the heat flux is applied. Sometimes a small wave 
develops just upstream of the heated zone. Ripples are apparent 
in the fluid over the patch heaters. Once the plate temperature 
exceeds a certain value, the film becomes very unstable and any 
small disturbance can trigger dryout. If no disturbance occurs, 
rupture is eventually induced purely by heat input. Dryout is very 
rapid when induced by heat input and spreads over large portions 
of the plate in a very short time. This mechanism is probably 
controlled by surface tension since this parameter is much larger 
in water than in the other fluids used. 

The second mechanism observed was dryout caused by 
instability due to surface waves. This mode is most apparent 
when using ethanol and Freon. Dryout occurs as follows: When 

heat is applied to the film waves appear on the surface. The film 
is thinner at the troughs of these waves so the evaporation rate is 
greater there. Vapor leaving the surface exerts a reactive force 
down onto the fluid which tends to deepen troughs and increase 
wave amplitude. Waves continue to grow until the film is thin 



enough in some area that a dry spot occurs. This dry area spreads 
until equilibrium is reached. 

The third mode is dryout due to gradual thinning. This 

mechanism was observed with ethanol, methanol and Freon. In 

this case, when heat flux is applied to the film, it begins to thin in a 
progressive manner. This thinning goes on until only a micro-thin 
layer of fluid is left on the plate. The thin micro-layer then 

spreads in the same way as a dryout. In fact, it is sometimes 

difficult to tell that the plate is still wet. If heat is applied long 
enough this micro-layer can also be burned off, however, this 
makes rewetting the plate much more difficult. This dryout mode 
seems to be dependent on inter-molecular forces such as London 
and Van der Waal's forces. 

Distilled water and tap water proved to be the most difficult 
fluids with which to create thin films. Water did not automatically 
wet the plate, as did the Freon TF. Water runs down the plate in 
small irregular streams unless forced to do otherwise. To create a 
thin film with water, it was necessary to first flood the plate and 
then slowly drain away excess fluid until only a thin layer 
remained. 

Once dryout occurred with water, it did not spread outward 
in a uniform way. Instead, the dried out regions moved about in a 
chaotic fashion. With standard tap water, dryout would usually 
encompass the entire heated zone, as would be expected. 
However, when distilled water was used as the working fluid, 



dryout would spread over the entire plate in a very short amount 
of time. Thin films formed with distilled water were found to be 
the most unstable. 

Methanol was another fluid which would not always wet the 
plate of its own accord. But here, the fluid could be made to form 
a thin film by wetting dry regions with a squirt bottle. Films of 
methanol were more stable than those formed with either tap or 
distilled water. They could only be ruptured through application 
of a heat flux. The other alcohol that was used was denatured 
ethanol. This fluid wet the plate naturally. To form a thin film of 
ethanol, it was only necessary to fill the reservoir at the leading 
edge of the plate. This was also the case for Freon TF, which was 
the fluid most willing to wet the plate. 

Those fluids which did not willingly wet the surface of the 
plate were also the ones which had the most trouble rewetting the 
plate after dryout had been caused by a heat input. Neither tap 
water nor distilled water would ever rewet the plate, even after 
the entire system had returned to room temperature. Methanol 
would rewet sometimes, but other times it would leave behind 
small patches that could not rewet naturally. Methanol was 
generally unpredictable in its behavior and was very sensitive to 
factors such as impurities on the plate or ambient humidity. For 
this reason, it was not used in acquiring data. Distilled water was 
not used because of its extremely unstable film characteristics. 


Freon TF and ethanol both rewet the dried areas in the same 
gradual fashion although Freon did so more rapidly. 

The contact angles exhibited by the various fluids varied 
greatly. Although no provision was made to accurately measure 
contact angles, differences were obvious through observation. 
Distilled water and tap water formed contact angles in excess of 90 
degrees while Freon TF maintained almost a zero contact angle. 
The contact angles for ethanol and methanol were also very small. 

Perhaps as a result of its large contact angle, water did not 
form micro-layers. Dryout was complete in all cases, with no 

residue being left behind on dried areas. On the other hand. Freon 
TF and ethanol did form micro-layers. This behavior is not the 
same as that observed by Mirzamoghadam and Catton [72] in their 
studies of meniscus behavior. 

In order to obtain a complete set of data, experiments were 
conducted while maintaining constant flow rate, constant film 
thickness, or constant angle of inclination for each of the fluids; 
water, denatured ethanol, methanol, and Freon TF. The 

thermophysical properties of these fluids are listed in Table 6.1. 
The method for maintaining the film at a constant thickness while 
varying the angle of inclination was as follows. At any given angle, 
increasing the level of head in the upper fluid reservoir produced 
a corresponding increase in the flow rate and film thickness. Thus, 
at each of the three angles of inclination used, the head in the 
upper reservoir was adjusted to produce the constant film 
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Table 6.1 Thermal Properties of the Working Fluids 


Fluid 

Freon TF 

Ethanol 

Methanol 

Water 

Viscosity 

(cp) 

0.682 

1.773 

0.82 

1.798 

Density 

(g/cm z ) 

1.56 

0.791 

0.810 

1.00 

Surface 

Tension 

(dynes/cm) 

17.3 

24.05 

24.49 

73.05 

Latent 

Heat 

(xltfj/kg) 

0.147 

0.949 

1.21 

2.68 

Thermal 

Conductivity 

(W/mK) 

0.070 

0.182 

0.215 

0.650 




thickness desired. The results for constant flow rates are shown in 
Table 6. II, while Tables 6. Ill and 6. IV show the results for 
constant film thicknesses and constant inclination angles, 
respectively. For accuracy and repeatability, all film thickness and 
flow rate measurements were made while maintaining the plate 
temperature at 25°C. 

Co ns ta nt Flo w R ato 

For the cases using a constant flow rate, the magnitude of the 
flow was different for each fluid. Because of variations in 
viscosity, wettability, and dryout heat fluxes, each fluid had a 
different flow rate which could be kept constant while varying the 
angle of inclination. Water required a flow rate of 0.89 ml/s while 
ethanol only required 0.26 ml/sec. A flow rate of 0.58 ml/s was 
required for Freon TF. While keeping the flow rate constant. 
Freon TF and ethanol exhibited similar patterns in that both the 
heat flux and plate temperature at dryout increased as the 
inclination angle decreased, as shown in Table 6. II. The opposite 

was true for water at constant flow rate which required an 
increase in heat flux and plate temperature as the angle was 
increased. This happened in spite of the fact that the film 
thickness decreased with increasing angle for all three fluids. The 
dryout heat fluxes for Freon TF ranged from 1.06 W/cm 2 on the 
horizontal plate to 0.45 W/cm 2 on the plate with a 5 degree 
inclination with plate temperatures at dryout ranging from 68°C to 
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Table 6. II Dryout and Rewet Heat Fluxes; Constant Flow Rate 


Fluid 

Flow Rate 
(ml/sec) 

Thickness 

(mm) 

Angle 

» 

dryout 

wycm 2 

T 

dryout 

-c 

n 

H rewet 

W/cm 2 

T 

rewet 

-c 

Freon TF 

0.58 

0.47 

0.0 • 

1.06 

68 

0.53 

57 



0.27 

2.5 ■ 

0.51 

50 

0.37 

48 



0.20 

5.0 • 

0.45 

45 

0.25 

41 

Ethanol 

0.26 

0.58 

0.0 - 

1.27 

78 

0.87 

68 



0.28 

2.5 * 

1. 13 

70 

0.81 

62 



0.24 

5.0 - 

1.04 

67 

0.71 

56 

Water 

0.89 

0.52 

0.0- 

0.95 

61 

- 

- 



0.31 

2.5- 

1.55 

81 

- 

- 



0.27 

5.0- 

1.83 

92 

- 

- 





























































4 5 °C. Similarly, ethanol’s dryout heat fluxes and plate 
temperatures ranged from 1.27 W/cm 2 to 1.04 W/cm 2 and 78°C to 
67 °C, respectively. The values for water, however, increased from 
0.95 W/cm 2 on the horizontal plate to 1.83 W/cm 2 on the plate 
inclined at 5 degrees and at 92°C. 

Rewet heat fluxes for ethanol and Freon TF followed a 
pattern similar to that of dryout. Increasing the angle of 
inclination led to a decrease in the heat flux which would allow 
rewet to occur. As with dryout, this seems reasonable since higher 
angles produce thinner films if the flow rate is held constant. With 
Freon heat fluxes for rewet ranged from 0.53 W/cm 2 at 57°C to 
0.25 W/cm 2 at 41°C. Ethanol rewet the horizontal plate at 0.87 
W/cm 2 at 68°C whereas 0.71 W/cm 2 at 56°C was the rewet limit 
on the 5 degree inclined plate. As mentioned before, water did not 
rewet the plate under any circumstances. Therefore, no rewet 
data is available in this case. 

Constant Film Thickness 

The flow rates needed to maintain constant film thickness, 
along with corresponding dryout and rewet heat flux data, are 
shown in Table 6.III. Thicknesses of approximately 0.3 mm were 
used for the different fluid films. All three fluids behaved 
similarly under the constant thickness constraint. As the angle of 
the plate was increased, the heat flux needed for dryout also 
increased. For Freon TF and water, the plate temperature also 
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Table 6. Ill Dryout and Rewet Heat Fluxes; Constant Film Thickness 


Fluid 

Flow Rate 

Thickness 

Angle 

ri 

*/ dryout 

T 

dryout 

» 

rewet 

T 

rewet 1 


(ml/sec) 

(nun) 


iV/cm 2 

•c 

W/cm 2 

•c 

Freon TF 

0.06 

0.29 

0.0 • 

0.32 

38 

0.19 

35 


0.71 


2.5" 

0.42 

45 

0.32 

43 


111 


5.0 * 

0.61 

49 

0.37 

44 

Ethanol 

- 

0.30 

0.0 • 

- 

- 

- 

- 


0.28 


2.5" 

1.22 

76 

0.91 

67 


0.36 


5.0 * 

1.29 

74 

0.91 

62 

Water 

- 

0.31 

0.0 • 

- 

- 

- 

- 


0.67 


1.2' 

0.57 

51 

- 

- 


0.93 


2.5" 

1.55 

80 

- 

- 





















































increased with increasing angle and heat flux, while for ethanol, 
the plate temperature for the higher heat flux case was actually 
less at dryout than for the lower heat flux case. This was also the 
case for rewet temperatures. In examining films of water and 
ethanol at constant thickness, the horizontal angle of inclination 
was not included. Even with near zero flow rates, the horizontal 
films could not be made as thin as the thickest film achievable at 
five degrees inclination. Thus one constant thickness could not be 
produced at all three plate angles being observed. For this reason, 
data for water and ethanol at constant film thickness includes only 
the 2.5 and 5.0 degree inclinations. Heat fluxes and plate 
temperatures for dryout varied from 0.32 W/cm 2 at 38°C for 
Freon TF on the horizontal plate to 1.55 W/cm 2 at 80°C for water 
on the 5 degree inclined plate. The lowest values for rewet were 
0.19 W/cm 2 at 35°C which occurred on the horizontal film of Freon 
TF, whereas the highest values of 0.91 W/cm 2 at 67°C occurred on 
a film of ethanol at a 5 degree inclination. 

Constant Angle of Inclination 

Table 6. IV gives the results for fluid films at constant angles. 
These results simulate a wicking material which exerts a constant 
capillary force on the working fluid. All three fluids exhibit the 
same pattern of increased heat fluxes and plate temperatures 
required for dryout as the flow rates and thicknesses are 
increased. Freon TF and ethanol were both run at 2.5 degree 
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Table 6.IV Dryout and Rewet Heat Fluxes; Constant Inclination Angle 


Fluid 

Flow Rate 
(ml/sec) 

Thickness 

(mm) 

Angle 

it 

dryout 

W/cm 2 

T 

dryout 

* c 

n lf 

* rewet 

W/crrr 

T 

rewet 

•c 

Freon TF 

0.43 

0.27 

2 . 5 • 

0.38 

45 

0.28 

43 


0.71 

0.28 


0.42 

45 

0.32 

43 


0.80 

0.30 


0.51 

50 

0.38 

47 

Ethanol 

0.31 

0.27 

2 . 5 ' 

0.71 

54 

0.48 

52 


0.52 

0.30 


0.91 

63 

0.68 

56 


0.71 

0.31 


1.08 

70 

0.77 

62 

Water 

0.74 

0.31 

1 . 2 ' 

0.57 

51 

- 

- 


0.93 

0.33 


0.75 

55 

- 

- 


1.79 

0.45 


1.06 

64 

- 

- 
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inclines whereas water was run at a 1.2 degree slope. The 

decreased angle for water was necessary so that a wide variety of 
flow rates could be examined while still allowing dryout, given our 
heat flux limit. The thicknesses of all of the films were 
approximately 0.3 mm with the exception of water at high flow 
rate, in which case the film was 0.45 mm thick. The lowest values 
of heat flux and plate temperature for dryout were 0.38 W/cm 2 at 
45°C for Freon while the highest values were 1.08 W/cm 2 at 70°C 
for ethanol. Rewet heat fluxes ranged from 0.28 W/cm 2 to 0.77 
W/cm 2 at 43°C and 62°C, respectively. 

6.5.3 Conclusions 

The various fluids used to form thin films displayed very 
different characteristics under dryout conditions. Although dry 
patches were formed through input heat flux, the stability of the 
dry patches was largely affected by the liquid contact angle. The 
problem is dynamic, since dry patches are constantly changing 
shape. One must therefore consider the moving contact line. A 
generally accepted formula for predicting dynamic contact angles 
has yet to be proposed although Orell and Bankoff [71] proposed 
two methods for approximating a static contact angle. 
Experimental measurements of dynamic contact angles usually 
involve direct observation through magnification. 

Some dimensionless parameters related to the dynamic 
contact angle were suggested by Dussan [73]. These include the 
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Reynolds, Bond, Weber, and capillary number, and the static 
contact angle. Contact angle hysteresis must also be considered. 
Hysteresis is affected by surface roughness and is not easily 
predicted. 

Wettability of a fluid on a solid surface is not determined 
solely by the contact angle. It was found that although ethanol 
easily rewet the copper plate used in the experiment, methanol 
would not naturally wet the plate. But the contact angles for these 
two fluids appear similar upon observation. Using the formula 
suggested by Orell and Bankoff [71] the static contact angles for 
the two fluids are within one percent. There seem to be many 
factors contributing to dryout and rewet, and none of them can be 
considered separately. 

6.6 CONCLUDING REMARKS 

Transient processes in the vapor core of heat pipes are very 
fast (few seconds). Therefore, the mentioned optical methods 
which provide inertialess measurements, are the proper tools to 
obtain accurate results. The temperature and concentration 
profiles are evaluated from the holograms developed by the twin 
beam interferometry. In this method, the temperature is evaluated 
from the difference between the phase shifts corresponding to two 
different wavelengths (see Eq. (6.11)). This difference is usually 
very small, therefor, the two wavelengths used should be as far 
apart as possible. 
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In laser Doppler anemometry, seeding of particles with 
proper density and size is a challenging problem. In systems like 
heat pipes, in which the working fluid undergoes a phase change, it 
is not possible to have particles recirculating with the fluid flow. In 
this case, particles can be carefully injected into the vapor flow in 
the evaporation end and collected from the condensation end. But 
the disadvantage of this technique is that the measured velocities 
at both 

ends are not accurate. 

When the evaporation and condensation rates are high (more 
than 1 w/cm 2 ), small liquid droplets with radius of order of 1 pm 
are present in the vapor flow. These droplets flow with the vapor 
and are able to scatter the laser beam. Therefore, no seeding is 
required in this case. 

There are several difficulties with the use of optical methods 
in heat pipe vapor flow measurements. Although these methods 
are very accurate and reliable, they are however, very expensive. 
The optical components are also very sensitive and not durable. 
Furthermore, the alignment of the laser beams passing through 
different parts of the measurement setup (specially in the laser 
Doppler anemometry), is hard to be achieved and easy to be lost. 

The major problem in the experiments, on the other hand, is 
the condensation of the vapor on the side walls of the rectangular 
pipe. A liquid film is formed on the walls, which disturbs the laser 
beam passing through the vapor core, and causes beam reflection. 


The unwanted condensation film was reduced by making the side 
walls out of acrylic plate which has very low thermal conductivity. 
The heat lost through the walls was highly reduced and, 
consequently, the condensation was confined to a smaller area 
closer to the liquid layer. In order to further reduce the 
condensation film, the Crystal Blue anti-fog jelly was coated on the 
internal surfaces of the acrylic walls. 



Chapter VII 

SUMMARY AND CONCLUSIONS 


The transient behavior of a heat pipe vapor flow is an 
important component of its operation and response to variable 
loading conditions. However, a simple model is needed in the 
thermal management of systems which use heat pipes as the heat 
transfer medium. A comprehensive experimental and theoretical 
research is needed to backup and verify the simple model. This 
study investigates the detailed dynamics of the vapor flow in heat 
pipes. A complete model will combine this work with studies of the 
liquid phase. 

The vapor flow in heat pipes is modeled and the governing 
equations are derived for transient two-dimensional compressible 
flow along with inflow and outflow boundary conditions. The 
equations are solved numerically for two cases; operating and 
startup transients. In the study of the steady state and operating 
transients after a sudden increase in the input heat flux, it is 
shown that for high-input heat fluxes, flow reversal occurs in the 
adiabatic and condensation regions. The circulation flow is caused 



by injection and extraction of a high momentum flow in the 
evaporator and condenser, respectively. 

During the startup transient in heat pipes, if the initial static 
state is at a very low pressure (which is most likely in aerospace 
applications), a high-input heat flux will create a shock wave in the 
evaporator. A nonlinear filtering technique is therefore used to 
capture the shock and overcome the cell Reynolds-number 
stability problem as a result of high velocity flow. Multiple shock 
reflections in the evaporator causes periodic oscillations in the 
pressure and, consequently, the pressure drops along the heat 
pipe. Depending on the transient time constant of the liquid layer 
and how fast the input heat reaches the vapor, the pressure drop 
in the vapor core may be significantly higher than its value in the 
steady state operation for the same input heat flux. Furthermore, 
the shock reflection causes a pressure increase in the evaporator. 
If the pressure exceeds the saturation pressure, mist nucleation 
will occur which will affect the evaporation rate. 

The optical methods are found to be the most appropriate 
and unique for temperature and velocity measurements in the 
vapor flow in heat pipes. Although these methods are very 
accurate, they are very sensitive and difficult to work with. The 
alternative approach would probably be to build a cylindrical heat 
pipe and to measure the input heat flux and the corresponding 
surface temperature after a sudden increase in the input heat flux. 
The standard and simple measurement technique, like 



thermocouples, could be used. Then the experimental results 
• would verify the results of the analysis which combine the vapor 

and liquid phases as well as the heat pipe case. This study is left 
for future efforts. 
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